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COMPUTER PROGRAM FOR THE TRANSIENT RESPONSE OF 
ABLATING AXISYMMETRIC BODIES INCLUDING THE 
EFFECTS OF SHAPE CHANGE 

By Lona M. Howser and Stephen S. Tompkins 
Langley Research Center 

SUMMARY 

A computer program to analyze the transient response of an ablating axisymmetric 
body including the effects of shape change is presented in detail. The program, its sub- 
routines, and their variables are listed and defined. The computer input and output, in 
printed and plotted form, for three sample problems are presented to aid the user in 
setting up and running a problem with the program. The governing differential equation, 
the boundary conditions for the analysis on which the computer program is based, and the 
method of solution of the resulting finite -difference equations are discussed. 

INTRODUCTION 

A numerical analysis of the transient response of an ablating axisymmetric body 
including the effects of shape change is presented in reference 1. The present paper 
briefly describes the analysis in reference 1 and presents in detail the associated com- 
puter program (program D2430) developed at the Langley Research Center. This paper 
also provides the user with an operating manual for the program. 

Some of the features of the analysis and the associated program are (1) the ablation 
material is considered to be orthotropic with temperature -dependent thermal properties; 
(2) the thermal response of the entire body is considered simultaneously; (3) the heat 
transfer and pressure distribution over the body are adjusted to the new geometry as 
ablation occurs; (4) the governing equations and several boundary -condition options are 
formulated in terms of generalized orthogonal coordinates for fixed points in a moving 
coordinate system; (5) the finite -difference equations are solved implicitly; and 
(6) other instantaneous body shapes can be displayed with a plotting routine. 

The computer program is written in the FORTRAN IV language for the Control 
Data 6000 series digital computer with the SCOPE 3.0 operating system. The equations 
have been programed so that either the International System of Units or the U.S. Custo- 
mary Units may be used. 



SYMBOLS 


A 


j_ 95 
x b 94 


A 


c 


VV c r l, i 

A s 


B s 

C 

C P 

H 

AH C 

AH S 

h l» h 2’ h 3 


K 

k 

L 

M 



m,n 


constant in oxidation equation corresponding to specific reaction rate 
coefficients in equations (6) 
constant in sublimation equation 

constant in exponential of oxidation equation corresponding to activation 
energy- 

constant in exponential of sublimation equation 

oxygen concentration by mass 

specific heat 

total enthalpy 

heat of combustion 

heat of sublimation 

coordinate scale factors (eqs. (2)) 

reaction-rate constant for oxidation (eq. (10)) 

thermal conductivity 

number of stations in x -direction 

molecular weight of gas 

molecular weight of oxygen 

integers 
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m 

“c 

“s 

p 

p w 

<*c 

^C,net 

q r 

R 

R cyl 

R stag 

r 

S 

T 

t 

w, z 

x, y 

x b 

a 


mass loss rate 

mass loss rate due to combustion 

mass loss rate due to sublimation 

exponent of pressure in sublimation equation (eq. (12)) 

wall pressure 

convective heating rate to nonablating cold wall 

hot -wall convective heating rate corrected for transpiration (eq. (9)) 

radiant heating rate 

radius of curvature of base curve 

cylindrical radius from axis of symmetry to base curve 
stagnation-point radius of curvature 

exponent of radius in sublimation equation (eq. (12)); spherical coordinate 

number of stations in y -direction 

temperature 

thickness of heat sink 

Cartesian coordinates (sketch 2) 

curvilinear coordinates (sketch 1) 

length of base curve 

absorptance 

weighting factor for transpiration effectiveness of mass loss due to 
combustion 
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o> weighting factor for transpiration effectiveness of mass loss due to 

sublimation 

(3 either 0 or 1 depending on whether transpiration or ablation theory is used 

6 material thickness 

e emittance 

|,?7 dimensionless curvilinear coordinates (eqs. (3)) 

9 angle between R and R ^ (sketch 1); spherical coordinate 

X mass of char removed per unit mass of oxygen 

p density of material 

a Stefan-Boltzmann constant 

r time 

p angle between axis of symmetry and normal to surface (sketch 1) 

Subscripts: 

e edge of boundary layer 

w wall condition 

x,y coordinates 

1,77 dimensionless coordinates 

Superscripts: 

* condition along x = L 

” condition along y = 0 
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DESCRIPTION OF MODEL 


Physical Model 

The analysis considers an axisymmetric ablating body exposed to aerodynamic 
heating; this body is composed of a single orthotropic material of varying thickness with 
temperature -dependent thermal properties. (See sketch 1.) The back surface of the body 
may be considered as a thin heat sink and/or radiator. Two coordinate systems are used 
to study the thermal and ablative response of the body. One is a curvilinear coordinate 
system, with x,y coordinates (sketch 1), which is used to determine internal temperature 
distributions. A stationary base curve located at the back surface of the body establishes 
the x-axis. 



The second coordinate system (sketch 2) is used to define the exterior geometry of 
the body which changes with time as a result of ablation. This coordinate system, with 
w,z coordinates, is a Cartesian system with the origin fixed at the original stagnation 
point of the body. All the geometric parameters needed to compute changes in the stag- 
nation heating rates and the heating -rate and pressure distributions over the surface are 
defined in this system. 

The governing time -dependent heat-conduction equation with variable coefficients 
for an axisymmetric body is, in fixed coordinates, 


1 

d l h 2 h 3 

h l h 2 h 3 

9x \ h i 



' h l h 3 


9 y\ hn 



= pc. 


9T 
P dT 


(1) 
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where the coordinate scale factors are 


h, = 1 + 

1 R 

(2a) 

h 2 =1 

(2b) 

h 3 = R cyl + y cos 9 

(2c) 


The transient temperature response of an ablating axisymmetric body is obtained 
from the solution of equation (1) with the appropriate boundary conditions, which are pre- 
sented in reference 1. The method of solution is discussed in the following section. 



Mathematical Model and Solution 

The finite -difference method was used to obtain the solution to equation (1). How- 
ever, if equation (1) were expressed in finite-difference form, it would describe the tem- 
perature variation at fixed stations in a fixed coordinate system. To maintain a fixed 
number of stations in a layer which changes thickness with time, it is necessary to change 
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the location of the stations and to interpolate to determine the temperatures at the new 
location after each time step. This procedure is time consuming and introduces a small 
error in each step of the calculation. This difficulty can be eliminated by transforming 
the equation to a coordinate system in which the stations remain fixed and the coordinates 
themselves move to accommodate changes in the surface location. 

This transformation can be made by introducing a moving coordinate system £, 77 , 

where 


and 77=1 
x b 6 


(3) 


In this system, the outer surface remains fixed at 77 = 1 and all other stations remain 
at fixed values of rj. 

The governing time -dependent heat -conduction equation (eq. (1)) in this transformed 
moving coordinate system is (eq. (9) in ref. 1): 

hhk 9 t\ 1 9 / h 3 . 9T\ _ J_ . nA 8T\ a f h 3 

8 2 37 ?^ 1 3 77 dv ) x fe 2 siyhj x b Z 5 dr]) 6 x b 977 ^ 


h l h 3 


9T' 


+ nA k 

6 2 5 877 \hj 977 


p \9t p6 dr] } 


(4) 


where 


A- 1 


(5) 


The unknown temperature field defined by the solution to equation (4) and its bound- 
ary condition was obtained by first approximating these equations by finite -difference 
equations with the use of the node pattern shown in sketch 3. Then the solution to these 
finite -difference equations is obtained with the method used in reference 2. 

This method is classed as an alternating-direction implicit method which has the 
advantages of being implicit, stable, and amenable to rapid solution. This method involves 
the alternate use of two finite -difference analogs to equation (1). In the first finite- 


difference equation at time t the analog to one of the second derivatives 


9^T 

9X 2 


, for 


example, is written at the new time r + A r, and the analog to the other derivative 


9 2 T 

9y 2 


is written at the old time r. Therefore, this equation is implicit in the x-direction (row) 
and explicit in the y-direction (column). 
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In the second finite-difference equation, at time t + 2 At, the analog 


written at the new time t + 2 At and the analog to 


9 2 T 

ax 2 


8 2 T 

ay 2 


is 


is written at the old time 


t + At. The second equation is implicit in the y-direction (column) and explicit in the 
x-direction (row). Using the two equations alternately results in a stable solution for any 
ratio of time increment to space increment as long as the same time increment is used 
for the successive application of the two equations. The time increment may be changed 
after the successive application of the equations. 


Equation (4) and the boundary conditions, when approximated by finite differences, 
lead to L sets of S simultaneous equations for a column solution and S sets of L 
simultaneous equations for a row solution. These equations take the form 


B 1 T 1 + C 1 T 2 ” D 1 


Vi-i + B i T j + c iVi = D i 


A N-1 T N-1 + b n t n 


= D. 


N 


(2 = j = (N - 1)) (6) 


where N is equal to S or L depending upon which finite -difference analog is applied. 
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Since the coefficients of equations (6) form a tridiagonal matrix, this set of simul- 
taneous equations can be quickly solved for temperatures. The method of solution based 
on the Gauss elimination method is discussed in reference 3. 

The coefficients of equations (6) are temperature dependent. Therefore, an itera- 
tion on these coefficients is made to obtain a temperature solution. 

OPERATION OF PROGRAM 

The physical problem to be modeled with the analysis is described by the FORTRAN 
input variables listed in a subsequent section. For example, the external body geometry 
is described in the w,z coordinates (sketch 2) which correspond to the input variables RS 
and ZS; material density corresponds to the input variable RO; and the stagnation cold- 
wall heating rate corresponds to the input variable QCTAB, which is a time -dependent 
array. Other input variables are required which control the solution, specify boundary 
conditions, and determine output from the program. These variables are listed in a 
subsequent section. 

This section describes the various boundary conditions that are available and a 
plotting routine that may be used with the oulput. The computation of the computing 
interval is also discussed. 


Boundary Conditions Along Front Surface 
An energy balance at the surface is 


/ H\ 

r 

H P / 

/h \ 

2 

2 


1 - (1 - 01 


- 0.084( — 

vcy 

S (“c“c + a s“ s ) 


+ aq r + m c AH c = k y |^ + m s AH s + aeT w 4 (7) 

where the terms on the left of the equality sign represent energy input to the surface and 
the terms on the right represent energy dissipation at the surface. The energy input may 
be any combination of convective heating, radiant heating, and the heat resulting from 
combustion. 

This energy input is accommodated by the heat conducted away from the surface 
and any combination of the heat radiated from the surface and the heat absorbed by sub- 
limation. The quantity of energy involved in each process is specified by the values 
assigned to the FORTRAN variables associated with that process. For example, the 


- j3(a» ni + a< ni 
c c s s 
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FORTRAN variables associated with the radiant heating rate q r are QRTAB, ALP HAT, 
and QRRAT, all of which define the radiant heating to the body with time. 

The pressure and the convective and radiant heating rates are functions of the body 
shape and also vary over the body surface. The changes in and q r at the stagna- 
tion point and the changes in pressure, q^, and q r around the body are computed within 
the program by setting IADJUST to a value greater than zero and specifying values for the 
variables defining the flow field and the body geometry. If IADJUST equals zero, then the 
variation of q^, q r , and the pressure over the body are tabulated as QRAT, QRRAT, 
and PRAT, respectively. 

Equation (7) shows that the mass loss due to combustion m c and mass loss due to 
sublimation m s affect the energy balance. This effect can be specified by either trans- 
piration theory (|3 = 0) or linear ablation theory (j3 = 1) . 

The rates of mass loss by both oxidation and sublimation are computed at each 
time step. However, only the larger of the two is used. 

The rate of mass loss by combustion may be specified by a half-order or a first- 
order oxidation equation. The input XORDER specifies which equation is used. The 
equation for a half -order oxidation reaction is (eq. (15) in ref. 1) 


M w( H e - H w) k2 P 


m = i<_ :~ w . \: ' e 

c 2 ] 1VU „ A 


w 


1 0 2 q C,ner 


M w(»e ' H w) k2 P 


“i2 


M 0 2 q C,net X 


+ 4K 2 C e p 
M 0, w 


(8) 


where 


H 


q C,net ” q C " ^)\ 1 " ^ ^ 

' e> 


H /r \2 

0.6 —(a m + cxem,,) - 0.0841 — J la 
q c v c c s sj \ q Q J \ 


„m„ + 

c c s s) 


/ * . x/V 

- e{a c m c + a s m s )l — 


( 8 ) 


and 


K = A c e 



( 10 ) 


The equation for a first-order oxidation reaction is (eq. (16) in ref. 1) 
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( 11 ) 


m c = rrsLS 

“02 KP w ( B e- H w) 

“» + “c,„et X 

The rate of mass loss by sublimation is (eq. (17) in ref. 1) 



( 12 ) 


Boundary Conditions Along Back Surface and Edge of Body 

Several boundary conditions may be specified along the surfaces at y = 0 and 
x = x^. These conditions are a constant -property heat sink, radiation from these surfaces 
to a surface at a specified temperature, or any combination of these. A heat sink along 
the back of the body is specified by the inputs CPDP, RODP, TDPRIME; along the edge of 
the body, by CPP, ROP, and TPRIME. Radiation from these surfaces is specified by the 
inputs EPSONPP, EPSONEP, and TBTAB. 

Output Plotting Routine 

The plotting routine for this program is convenient for studying the results of calcu- 
lations. This routine is activated by setting IPLOT equal to an integer greater than zero. 
The following plots are generated: (1) RSS versus ZS at times listed in the PLTIME 
table (this plot shows the body geometry), (2) MDOT versus X at each PRFREQ time, and 
(3) T(N) versus X at each PRFREQ time, where N is a specified row of temperatures. 

For example, to plot the temperatures of rows 2, 6, and 8, set the input NTP = 3, 2, 6, 8, 
where the 3 specifies the number of rows to be plotted. Other input quantities that must 
be specified are MDMAX, RSSMAX, ZSMAX, PTMAX, and PTMIN. These inputs specify 
maximum and minimum values which are used to get reasonable plotting scales. Sample 
plots are shown with example problems discussed in a subsequent section. 

The plotting routines used are from the CalComp software package. Plotter output 
is routed to a tape during job execution and after job completion is plotted on a CalComp 
digital incremental plotter. 


Computing Interval 

Although the alternating-direction implicit method used for solution of the finite - 
difference equations has the advantage of being stable for any time increment, the choice 
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of a computing interval is important. An initial and a maximum computing interval 
DSLTAU and DTMAX are inputs for the program. After the application of a column and 
a row solution, the program computes an interval for the next two successive time steps. 
This is done by examining the number of iterations necessary for convergence at the pre- 
vious time step. If this number was (a) equal to 1, the computing interval will be doubled, 
but will not exceed DTMAX; (b) equal to 2, the interval will not be changed; or (c) equal 
to 3, the interval will be halved. 

This should not be confused with the input MAXITT. If the number of iterations 
during a solution that is not a row solution exceeds MAXITT, the computing interval will 
be halved and the solution restarted. 

PROGRAM DESCRIPTION 

The computer program D2430 was written in FORTRAN IV language for the Control 
Data 6000 series digital computer under the SCOPE 3.0 operating system. The program 
requires approximately 70 000 octal locations of core storage. 

This section presents the program, its subroutines, and their variables. The vari- 
ables are grouped in labeled COMMON blocks PICK, INPUTS, and HOLD. Input data are 
loaded with FORTRAN IV NAMELIST. The variables in INPUTS (except the variable 
DUMMY) and in HOLD are also in the NAMELIST statement which appears in another 
section. 


Labeled COMMON 

r 

The following list contains the FORTRAN variables appearing in labeled COMMON 
and the dimensions of the array for each variable. The notation is in the form A(m,n). 


COMMON FORTRAN 
label variable 

PICK 

A(10,20) 

AA(20) 

AB(10,20) 

ALPHA(20) 

B(20) 


Description 


Elements in coefficient matrix for the column solution 
86 
9x 

Elements in coefficient matrix for the row solution 


a 

Major diagonal elements in coefficient matrix 
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COMMON 

label 

FORTRAN 

variable 

Description 

PICK 

BS1(10,20) 

Major diagonal elements in coefficient matrix for the 



column solution minus — term 

9t 


BS1B(10,20) 

Major diagonal elements in coefficient matrix for the 

AT 1 

row solution minus — term 

9 r 


C(10,20) 

Elements in coefficient matrix for the column solution 


CB(10,20) 

Elements in coefficient matrix for the row solution 


CK(10) 

Temporary storage used to define the thermal 
conductivity at a half station 


CKETA(10,20) 

at the station 


CKXI(10,20) 

at the station 


COST (20) 

cos 9 


CP(10,20) 

c p 


D(10,20) 

h 2 h 3 k | 

h l 


DC (20) 

Right-hand side of the matrix solution 


DELESQ 

(At]) 2 


DELETA 

At] 


DELTA(20) 

5 


DELXI 

AS 


DELXISQ 

(A£) 2 


E(10,20) 

h l h 3 k 7/ 

h 2 


EIGHT3 

Constant, 8. 0/3.0 


ELAM(20) 

X 


ETA(10) 

V 
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COMMON 

label 

FORTRAN 

variable 

Description 

PICK 

EXPG 

Computed constant used in computing new heating 



distribution 


F(10,20) 

h 2 h 3 k | ?? 95 
hj6 9x 


GG 

Computed constant used in computing new heating 



distribution 


GIMACH 

Computed constant used in computing new pressure 



distribution 


Hl(10,20) 

h l 


H2(10,20) 

h 2 


H3(10,20) 

h 3 


HC(20) 

AHg 


HCOMB(20) 

AH C 


HE 

H e 


HW(20) 

H w 


IFIRST 

Internal code; 0 for first time step in calculation, 1 for 



any time after first time step 


IROCOL 

Internal code; 1 for column solution, 2 for row solution 


ITC 

Number of iterations during the column solution 


ITR 

Number of iterations during the row solution 


ITT 

Number of iterations during a solution 


ITTO 

Total number of iterations from the initial time 


LM1 

Computed constant (L-l) 


LM2 

Computed constant (L-2) 


MCDOT(20) 

m c 


MDOT(22) 

m 


MSDOT(20) 

“s 
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COMMON 

label 

FORTRAN 

variable 

Description 

PICK 

PED2 

Constant 1.5707963268 


PRELOC(20) 

Local wall pressure 


QC(20) 

Adjusted convective heating rate 


QC1 

q C 


QCNET 

q C,net 


QCOMB(20) 

Heat due to combustion for oxidation 


QR(20) 

Adjusted radiant heating rate 


QR1 

q r 


QS(20) 

Net heat input 


RNS 

Nose radius 


RODPC 

t"p”c p ’7AT 


ROPCPP 

t'p'CpyAr 


RSS(22) 

Coordinate used to define body geometry, w 


RST02 

Computed constant, ratio of molecular weight of free 



stream to molecular weight of diatomic oxygen used 
in oxidation equation 


SIG 

Computed constant ae 


SIGDP 

Computed constant ae" 


SIGMA 

a 


SIGP 

Computed constant ae' 


SINT (20) 

sin 9 


SMI 

Computed constant (S-l) 


SM2 

Computed constant (S-2) 


TAU 

Time at which calculation is being made 


TB 

Temperature to which back surfaces radiate 


TT(10,20) 

Estimated temperatures at r 


TWDELXI 

Computed constant 2 A£ 
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COMMON 

label 

FORTRAN 

variable 

Description 

PICK 

TWOGI 

Computed constant used in computing new heating 
distribution 


V(20) 

Elements in coefficient matrix for column solution 


VB(10) 

Elements in coefficient matrix for row solution 


X(22) 

Curvilinear coordinate 


XDXISQ 

Computed constant (x^ 


XODXI 

Computed constant x^ A£ 


Y(10,20) 

Curvilinear coordinate 


Z(20) 

Elements in coefficient matrix for the column solution 


ZB(10) 

Elements in coefficient matrix for the row solution 

INPUTS 

DUMMY 

Used in setting initial values of all inputs to zero 


All the variables in NAMELIST except TMIN also appear in INPUTS 
HOLD TMIN A minimum temperature value 

Descriptions, Flow Charts, and Listings 

This section identifies the main program and each subroutine in the program D2430. 
A brief discussion, a flow chart, and a listing for each are given. The numbers appearing 
in the flow charts represent a FORTRAN statement number in the program. The interpo- 
lation subroutines FTLUP and DISCOT are described in detail in appendix A. 

Program D2430 .- Program D2430 is the control program. It reads the inputs, calls 
the subroutines to solve for the temperature profile, calls subroutines for plotting, and 
controls the iteration scheme for the temperature solution. The flow chart for pro- 
gram D2430 is given on the following pages: 
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The listing for program D2430 is as follows: 

PROGRAM 02430 ( I NPUT, OUTPUT , T APE5 = I NPUT , TA PE5=CUTPUT ,T AP £7=20 1 , 

1 TAPES =231, TAP £<5=201 I 

AX I S YMMETR I C ABLATION PROGRAM 

TWO-DIMENSIONAL ABLATICN ANALYSIS FOR AXIALLY SYMMETRIC BODIES OF REVOLUTION 
AT HIGH HEATING RATES j CONSIDERING SHAPE CHANGE 

THIS IS THE MAIN PROGRAM - IT CONTROLS THE GENERAL FLOW CF PROGRAM 

COMMON /PICK/ A( 1C, 20) ,AA(20»,AB( 10,20) ,ALPHA(20) ,B(20) , 

2 BS1U0.2C) ,BS18Q0,20),C(i0,20),CB(lC,20) ,CM10 > ,CKETA( 10,20 > , 

4 CKXI( 10,20) , COST I 20) ,CP ( 10 , 20 I ,D ( 1C , 20 ) , DC < 20 I , 

6 DELESQ,0ELETA,DELTA<20) ,CELX I , DELX I SC, E < 10 ,20 ) .EIGHT3, 

8 EL AMI 20) . ETA ( 1C ) ,EXPG ,F (1 0 , 20 I , GG, GIM ACH, HI ( 10 , 20 ) , H2 ( 1 C, 20 I , 

A H3I10 ,20) , HC ( 20 ) , HC0M8I 20 ) , HE , HW ( 20 I , IF IR ST, IROCOL , I TC , ITR , I TT, 

C ITT0,LM1 ,LM2, MCDOT ( 20 ) , MC OT ( 22 I ,MSD0T(20) , P 102 , PPELOC ( 20 ) , QC ( 20 ) , 

E CC1 , QCNET (20 ( ,CCCMB ( 20 ) ,QR ( 20 ) , QR 1 , QS( 2C I , RNS , RODPC ,ROPCPP , 

G R SSI 22) ,RST02,SIG,SIGDP, SIGMA, SIGP,SINT(20 ) , SMI , SM2 ,T AU ,TB, 

1 TTI10 ,20 ) ,TTF ( 10 , 20 ) , TWQE LX I , TWOGI ,V(20),VB(10),X(22),XDXISQ, 

K XCDX I , Y ( i 0 ,20 ) , Z 1 20 ) , ZB 1 3 0 ) 

COMMON /INPUTS/ DUMMY! i) , AEXP , ALC TAB U 0 ) ,TTALC(10) , M AL PHC, NAL PHC , 

2 ALPHATI10I ,TALPHA(10) ,M ALPHA , NAL PHA , ALST AB ( 10 I , TT AL SI 1C) ,MALPHS, 

4 NALPHS,ASEXP,BETA,BEXP,BSEXP,CE,CKETATBI50) ,TTCKETA(10) , 

6 ETATABI5) ,NCKETA,NETA,CKXITAB( 50 I ,TTCK X 1 1 10 ) , XI TAB I 5 ) ,NCKXI , 

8 NXI ,C0RDSY»CPDP,CPP,CPTAB!10) , TT ABCP 1 10 ) , MCP , NCP , DELT AO 1 20 I , 

A DELTAU,0FLTMIN,DTMAX,ELAMTB(28), TTELAM17) , PEL AMI 4 ) , NE LA M , 

C NPELAM,ENDTAU,EPSCNE,EPSCNEP,EPSCNPP,ERRORT,GAMBAR,GAMINF, 

E HC0M3TB (28) , TTHCOMB I 7 ) , PHCOMB I 4) , NHCGMB .NPHCCMB ,HCT AB 1 2 8) , 

G TTABHCI7I ,PHC(4) ,NHC , NPHC , HET A8I 10 ) , TT ABHE I 10 ) ,MFE , NHE , HWTABI 15), 

1 TTABHWI 151 ,MHW,NHW, I A 0 JUST , I P LOT , L , MACHNO , M AX I TT , MDMAX , 

J MDCT0I20), 

K MW02, MWSTR,NTP(7) .PLTIME 1 15 ), PRAT 1 20) , PRFREQ , PSEXP , PSTAGTB 1 1 0 ) , 

M TT PST AG DC I , MPST AG, NP ST AG , PTMAX , PTM IN ,QCTAB I 1C) , T TABQC 1 1C ) ,MQC, 

N NCC , QRAT I 20 ) , 

0 GRRAT (20) .QRTABUO) .TTABORI 10 ) ,MQR ,NQR ,R! 20 ) ,RI EXP ,RN SI ,RO,RCDP, 

Q RGP,RS(2CI ,RSSMAX,S,STEBCL,T( 13,201 , T AUC , TBT AB (10 ) , TT ABTBI 10 ) , 

5 MTB,NTB,TDPPIME,THETA(20) ,T PR IME , XO, XCRDE R, ZS ( 20 ) , ZSMAX 
DIMENSION DELTdO ,20) ,ZZ(22I ,Y3L(2) 

REAL MDOTC , MDO T , MCOOT , MSDGT , MW STR , MW02 , MACHNC, MDMAX 
INTEGER S , SMI , SM2 

DATA XLABFL,YLABEL,X2L ,Y2L,Y3L/ 2HZB , 3HR S S , 1 HX ,4HMD0T , l 2HTEMPERATU 
IRES/ 

NAMELIST /D2430/ AEXP , ALCTAB , TT ALC ,M ALPHC , NAL PFC , A LPHAT , 

2 TALPHA.MALPHA , NAL PHA, ALST AB , T T AL S , M AL PHS , NAL PHS , A SE XP , 

4 BE TA, BE XP ,B SEXP ,CE ,CKET ATB , ET ATA B ,TTCK ETA, NCKETA.NET A, 

6 CKXITA8,XITAB,TTCKXI,NCKXI, NXI, CORDSY,CPDP,CPP,CPTAB,TT ABCP, MCP, 

8 NCP.DELTAQ, DELT AU ,DELTM IN , OTM AX , ELAMTB , TTEL AM ,PEL AM ,NEL AM, NPELAM , 

9 ENOTAU.EFSCNE .EPSONEP ,E PSCNPP .ERROR T, GAMBAR ,G AMI NF ,HC CM BTB , 

A TTHCOMB, PHCOMB, NHCOMB ,N PHCOMB, HC TAB , TT A8HC , P HC , NHC , N'PHC ,HE TAB , 

C TTABHE,MHE,NHE,HWTAB,TTABHW,MHW,NHW,I ADJUST, I PLOT, L, MACHNO, 

E MAXITT,MCMAX,MD0T0,MW02,MWSTP ,NTP , PLT IME , PRA T , PRFREC , PSEXP , 

G PSTAGTB, TTPSTAG.MPSTAG.NP STAG, PTMAX,PTMIN,QCTAB,TTABQC, MQC, 

1 NQC,QRAT,QRRAT,QRTAB,TTABQR,MGR,NQR,R,RIEXP,RNSI,RO,RCOP,ROP,RS, 

K RSSMAX, S, STEROL, T, T AUO.TB TAB, TTABTB,MTB,NTB,TDPRIME, THETA, 

M TMIN,TPRI.ME,XC,XORDER,ZS, ZSMAX 
CCMMCN /HCLD/ TMIN 
TM IN =0. 

DC 10 1=1, S34 
10 DUMMY ( I )=C.C 
DTMAX=2, 

1 READ 15, ICC) 

100 FORMAT ( 80 F 

1 ) 

IF (EOF, 5) 2,3 

2 STOP 

3 READ ( 5 , D2430 ) 

WRITE ( 5,C2430 ) 

WRITE (6, ICC) 

C 

C SET INITIAL VALUES 


1C0000 
200000 
3C0000 
4 COOCO 
500000 
bOOOOO 
700000 
8C00C0 
900000 
ICCOnoO 
1100000 
1200000 
13000C1 
1400000 
1500000 
1600000 
1700C00 
1800000 
190CGC0 
2CC0000 
2100000 
22C0000 
2300000 
2400000 
2600001 
2600002 
2700000 
2800000 
2 90'Q000 
3000001 
3300000 
32COOGO 
3300000 
3400000 
350C000 
36000C0 
3700000 
38C0000 
3900000 
4000000 
4100000 
42COOOO 
4300000 
44COOOO 
45C0GC0 
4600000 
4 7000 CO 
4800000 
4900000 
5000001 
53 COOOO 
5200000 
5300000 
5 4000 CO 
5500000 
56C0000 
5700000 
5800001 
5900000 
60Q0UC0 
61C0000 
6200000 
63C0000 
6400000 
8500000 
bfcOOOCO 
7000000 
7100000 
72C0CC0 
73COOCO 
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NNTP = NTP < 1 5 

PID2 = 1® 57C7963268 

TWO GI = 2s 0 / ( { GAM INF - 1.0! * MACHNO **2 ) 

EXPG = (GAM'BAR - 1.0)/ GAMBAR 
GIMACH= 1. / { GAMINF * MACHNO **2 ) 

GG= SQRTI EXPG * (1*0 + TWOGI) * (1.0” GIMACH)) 

GG= SORT (GG) * 2.0 

INCP=0 

I R0W=0 

ID T=1 

DT AUO= 1.0 

DTAUi=DELTAU 

IROCOL =1 

C WILL PRINT ONLY AFTER A COL. AND RCW COMPUTATION HAS BEEN MADE 
T AUOO= TAUO+ PRFREQ 
ITT0=3 
DO 11 M=1,S 
DC 11 N=1,L 
CELT ( M ( N ) = 1000. 

11 TT(M,N)= T(M,N1 
DELTAU=DELTAL'/2. 

TAU=TAUO+CELTAU 
IF IRST = 0 

ITT=1 

LM1 = L- i 
ALM1 = L Ml 
LM2 = L- 2 
SMI = S- 1 
SM2 = S- 2 
DELXI =1. /ALM1 
DELX = XO/ ALM1 
RST02 = MWSTR/MW02 
X ( 1 » =0. 

00 12 N=2 »L 

12 X ( N ) = X(N-l) + OELX 
DELETA = 1./SM1 
OELX ISO = DELX I **2 
OELE SQ= DELETA **2 
TWDELXI = 2.0* DELXI 
EIGHT3=8. C/3.0 

DO 18 M=1,S 
AM=M— 1 

18 ETA(M) =DELETA*AM 
SIGMA=STEECL 
SIG a SIGMA* EPSONE 
SIGP = SIGMA * EPSONEP 
SIGDP= SIGMA * EPSCNPP 
XOOXI = XO * DELXI 

RODPC = TDPRI ME*RODP * CPDP / DELTAU 

ROPCPP = TPRIME * ROP * CPP/ DELTAU 

RODT= RO/DELTAU 

XDXISQ = XC**2 * DELXI SQ 

DC 22 N=1,L 

MOOT ( N ) =MCOTO( N ) 

MCDCT(NI=MDCTO(N) 

MS DOT ( N ) =MDOTO ( N ) 

20 DELTA ( N ) = OELTAC(N). 

THETA! N)=. C 1745 32925*THETA INI 

SINT(N) = S IN ( THETAI N ) 1 

ZZ ( N ) = IS (M + CELTAO(N)*SINT(N) 

22 COST ( N ) = CCS ( TFtTA (N ) ) 

IF ( IPLOT, EC. 0 ) GO TO 23 
C PLOT BASE CURVE IF PLOTTING IS CALLED FOR 
REWIND 7 
REWIND 3 
REWINO 9 


7400000 

7500000 

7600000 

7700000 

7800000 

7900000 

8000000 

81C0000 

8200000 

83000C0 

8400000 

8500000 

8600000 

8700000 

88C00C0 

89COOOC 

9CC0C00 

91000C0 

9200000 

9300000 

9400000 

9500000 

9600000 

97COOOO 

98000C0 

9900000 

10000000 

10100000 

10200000 

10300000 

10400000 

10500000 

10600000 

10700000 

10800000 

10900000 

uccoooo 
11100000 
11200000 
11300000 
11400000 
11 5000C0 
11600000 
11700000 
118300C0 
11900000 
12000000 
12100000 
1 2200000 
12300000 
12400000 
12500000 
12600000 
127CC000 
128000C0 
12900000 
13CC0C00 
13300000 
13200000 
1330000C 
1 34C0000 
13500000 
13600000 
1370000C 
13800000 
139C0000 
14030000 
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non non on o 


CALL CALCCMP 
I PLT = 1 
I PL TK = 0 

IF ( CORDS Y® NE » C I GO TO 2250 
WRITE (7) (ZZ(N) »RS(N) ,N=1 ,L) 

GO TO 23 

2250 WRITE (71 ( ZS ( N ) ,DELT A <N ) , N=1 , L I 

COMPUTE H-S 

23 DC 25 M=1,S 
DO 25 N=1,L 

Y(M,N)=ETA(MI*GELTA(N) 

HI ( M * N I = 1.0 + ETA(M)* DELTA(M/R(N) 

H2 (M,N)=1. 

25 H3(M,N)= R S ( N ) + Y(M,NI *COST(N) 

05 DO 101 M = 1,S 
DO 101 N=1,L 

CALL FTLUP ( TT (M ,N I , CP ( M ,N I , MCP, NC P , TTABCP, CPTAB I 
CALL DISCCT (TT(M,N),X (N).TTCKXI , CKX I T AB , X I T AB , 1 1 , NCK X I , NX I , 
1CKXI (M,N) ) 

101 CALL DISCCT (TT ( M , N I , Y ( M,N) , TTCKET A , CK ET AT B , E T AT AB » 1 3 , NC KF TA » NETA , 
2CKET A( M,N I I 
AA ( 1 1 =0. 0 
DO 109 N=2»LM1 

109 AA(N)= (Of LTA( N+l) -DELTAIN-l ) ) /(TWDELX I*XO ) 

AA(L ) = ( 3. CODELTA (LI-4. 0*DE LTA( LM1 ) +DEL TA ( L M2 ) I/ITWDELX I*XO> 

DC 110 N=3 ,L 
DO 110 M = 1,S 

D ( M » N I = F2 (M,N)*H3(M,N)* CKX I ( M , N I /HI (M ,N I 

E(M,N>= HKK.N)* H3(M,N) * CKETA(M.N) / H2(M,N) 

110 F(M,NI=D(M,N)*fcTA(M)*AA(N)/OELTA(Nl 
CALL SQAERC 

GO TO (310,320 1 , IROCOL 
310 CALL COLUMN 
ITC= I T T 
IF IRST = 1 
GO TO 350 
320 CALL ROW 
ITR= ITT 

IF (IROW.EC.OI I R0W=2 
350 CONTINUE 

IF ANY TEMPERATURES ARE NEGATIVE STOP CALCULATIONS 
DO 360 N=3 ,L 
DO 360 M=1,S 

IF (TTF(M,M.LE,OI GO TO 411 
350 CONTINUE 

TEST TO SEE IF TEMPERATURES HAVE CONVERGED 

ITTO= I TTO+1 
DO 400 N=3 ,L 
DO 400 M= 1 , S 
ABSTT=ABS(7T(M ,NI I 
ABSTTF=ABS(TTF(M,NI) 

TEST=ABS( ABSTTF-ABSTTI/ABSTT 
IF (TEST - ERRCRTI 4CC,40C,7 00 
430 CONTINUE 

COMPUTE MDOT 

CALL SQ4ER0M 
C COMPUTE DELTA 

DO 410 N=i,L 

DELTA(N)=CELTAC(NI-(MCOTG(NI+MDOT(N) l*DELTAU/(2»0*RO) 

C RESET DELTAO AND MDCTO 
410 MDOTO( NI = MCOT(M 

C IF DELTA BECOMES LESS THAN DELTMIN (SOME MINIMUM DELTA INPUT) STOP 


14 1 OOC 00 
14200000 
14300000 
14400000 
14500000 
14600000 
147C000C 
14800000 
149C0000 
15C00Q00 
J.51C0000 
152COOOO 
15300000 
154QC000 
15500000 
15600000 
157C0000 
15800000 
15900C00 
16C0Q0C0 
16100000 
16200000 
16300000 
16400000 
16500000 
16600000 
16700000 
16800000 
169C0000 
17000000 
17100000 
17200000 
17300000 
17400000 
17500000 
176000C0 
177000CC 
17800000 
17900000 
18000000 
18100000 
18200000 
1 830000C 
184C0CCC 
185000C0 
18600000 
18700000 
18 8000C0 
18900000 
19000000 
19100000 
19200000 
193C00C0 
194C000C 
19500000 
19600000 
19700000 
19800000 
19900000 
20000000 
201000C0 
20200000 
20300000 
20400C0C 
205000C0 
20600000 
20700000 
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C THE C ALCUL AT I CNS 
DO 412 N = I,L 

IF ( DELTA ( N ) «GT . DELTM IN } GO TO 412 

411 CALL Z PR I NT 
STOP 

412 CONTINUE- 

IF ( INOP.EQ.l ) GO TO 418 
IF (TAU.LT.TAUCC) GO TO 420 
IF ( IROCOL. EQ. II GO TO 418 
IN0P=1 
GO TO 420 

418 I NOP =0 
TAUOO=TAUCG+ PRFREQ 

C 

C 

CALL ZPRINT 
C 

IF ( IPLOT . EQ. 0 ) GC TO 420 

I PLTK= IPLTK + 1 

WRITEI8) ( MDOT(N) , N=1 ,L) 

IF ( NNTP. EC. 01 GO TO 420 
DO 419 M= 1 1 NNTP 
1= NTPlM+1) 

419 WRITE (9) (TTF( I ,N) ,N=1 ,L1 

420 IF (IROW-II 540,490,484 
434 DELTAU =QELTAU*2,0 

IROW= 1 
KFRE=KFRE+1 
C 

C OBTAIN DEL T AU AS A FUNCTION OF ITERATION OF PREVIOUS TIME STEP 
490 DTAU1 = DFLTAIJ 

IF (IROCOL.EQ. 1) GQ TO 54G 
IF IITT-2) 49 5 , 54C ,530 
495 DELTAU=2.C*0TAU1 

IF (DELTAU.GT. DTMAXI DELTAU=DTMAX 
GO TO 540 

530 OELTAU=OTAll/2. 

IF (DELTAU.LT, l.E-6) GO TO 900 
540 TAUO = TAU 

C CHECK TO SEE IF IT IS TIME TC PLOT 
IF ( IPLOTr EC. 0) GO TO 543 
IF (TAU.LT.PLTIME( IPLT I) GO TO 543 
I PLT= I PLT+1 

IF ( CORDS Y . NE. 0 ) GO TO 542 
WRITE (71 (ZS(N) ,RSS(N),N=l,U 
GO TO 543 

542 WRITE (71 ( ZS ( M .DELTA (N », N*1 , L » 

C 

C INCREMENT TIME AND REPEAT CYCLE ALTERNATING ROW AND CCLUMN SOLUTION 

543 TAL-TAU+DELTAU 

RODPC = TOPRIME*RODP * CPDP / DELTAU 
RCPCPP = TPRIMF * ROP * CPP/ DELTAU 
ROOT = RO/CELTAL 
IF (TAU. GT.ENDTAU) GO TC 950 
C 

C EXTRAPOLATE TC GET NEW GUESS TEMP(TT) 

C 

00 446 M=i , S 
DO 445 N-l ,L 
DELT(.M,N» = 1000. 

DELTN=TTF(M,N)-T(M,N) 

T(M,N)=TTF(M,N) 

446 TT(M,N)=TTF(M,N) + (DELTAU/DTAU1 1 40ELTN 
GO TO ( 35C ,650 t , IROCOL 
550 IROCOL = 2 
ITT=1 
GO TO 25 


20800000 

20900000 

21000000 

21100000 

21200000 

2130000C 

21400000 

21500000 

21600000 

21700000 

21800000 

21900000 

22000000 

22100000 

222C0000 

22300000 

224C00C0 

22500000 

22600000 

227C0000 

22800000 

22900000 

23000000 

231COCOO 

23200000 

23300000 

23400000 

23500000 

23600000 

23700000 

23800000 

23900000 

24000000 

241C0000 

24200000 

24300000 

24400000 

24500000 

24600000 

24700000 

248000C0 

24900000 

25000000 

251000C0 

25200000 

253COOOO 

254C0000 

25500000 

256C00C0 

257C00C0 

25800000 

259COOCO 

26000000 

26100000 

262COOOO 

26300000 

26400000 

26500000 

26600000 

26700000 

268CC0C0 

26900000 

270C0000 

27100000 

27200C00 

27300000 

27400000 
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650 IROCOL = 1 
I T T = I 
GO TO 23 
C 

C TEMP. DOES MOT MEET ERROR CRITERIA, MUST ITERATE AGAIN 
C MEW GUESS IS TEMP, CF PREVIOUS ITERATION TT =TTF 
C 

700 ITT = I TT +1 

IF (ITT - MAXITT) 705,705,800 

705 00 720 M=1,L 
DC 720 M = 1 ,S 

DELT1 = ABS(TTF(M,M- TT(M,NI> 

IF (DELTl.LT.iO. ) GO TO 718 
IF (DSLT1 -CELT ( M , N) ) 718,750,750 

713 DELT(M,N)=CELT1 

7?0 CONTINUE 

DC 730 M=!,S 
DO 730 N=1,L 

730 TT ( M , N ) = TTF(M.N) 

GO TO 95 

750 IF (ITT.LT.31 GO TO 718 
C 

C PROGRAMED STOPS 
C 

WRITE (6,752) 

752 FORMAT (*CTEMPERATURE IS DIVERGING WHY*) 

758 WRITE (6,755) 

759 FORMAT ( *OTT ( M , N ) * I 
DC 765 M= 1 , S 
MM=S-( M-l ) 

765 WRITE ( 6,766 ) ET A( MM I , (TT( MM ,N ) ,N=1 ,L> 

766 FORMAT ( F6 ,3 , 6X1 5F8. 1/ ( 12X , 1 5F 8, 1 ) I 
WRITE (6,767) IROCOL 

757 FORMAT ( *f IROCCL = * I 3 ) 

CALL ZPRINT 
STOP 

300 IF ( IROCOL. EQ, 1 ) GO TO 8C3 
WRITE (6, 8C1) 

801 FORMAT (*OTHIS IS A RCW SOLUTION, DELTAU CANNOT CHANGE) 

GO TO 758 
C 

803 DT AU 1= DELTAU 

DELTAU = CELT AU/ 2« 0 
WRITE <6,fC5) DELTAU ,TAU 

305 FORMAT (*CI DID IT— DE LT AU=*E 14, 5 , *TAL=*E1 A. 5 I 

IF (DELTAU. LT. l.E-6) GO TO 930 
TAU = TAU - DELTAU 
DC 810 M=1,S 
DO 810 N=1,L 
DELT(M,Nt=100G» 

310 TT ( M ,M ) = T ( M , N ) 

ITT = 1 
GO TO 95 

900 WRITE (6,901 ) 

901 FORMAT ( +OTEMPERATURE ITERATION DOES NOT CONVERGE* ) 

GO TC 753 

C 

C PLOT ZS VS. RSS , X VS MOOT , X VS BACK SURFACE TEMPERATURE 
C 


27500000 
276C00CC 
277C.OOCO 
278G0000 
275D0000 
28G00C00 
28100CCD 
282Q0CC0 
28 330CC0 
28400003 
28500000 
28600000 
2870000 0 
288C00C0 
2890000C 
29CCOOCO 
2910000C 
29200000 
293000C0 
294CC000 
29500000 
29600000 
297000C0 
29800000 
2990000C 
300C0C0C 
30100000 
30200000 
3C300000 
304300C0 
30530000 
30600000 
307000CC 
30800000 
30900000 
31COOOCO 
31100000 
312009GC 
31 3D0CG0 
31400000 
31500000 
. 31600000 
31700000 
31 800000 
31900000 
32000000 
32100000 
322C0C00 
32300000 
324C0000 
32500000 
326CCOG0 
327C0000 
32800000 
32900000 
33000000 
33100000 
33200000 
333COOCO 
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959 CALL ZPRINT 

IF ( I PLOT* EC« 9 ) GO TQ 1 

END FILE 7 

END FILE 8 

END FILE 9 

REWIND 7 

REWIND 8 

REWIND 9 

IEC = 3 

DO 960 M= 1 » I PLT 

READ (7) IZZ(N), RSS(N) ,N=1 ,L) 

IF (M.EQ.IPLT ) IEC =1 

960 CALL INFOFLT ( I EC , L , Z Z , 1 , R SS , 1 , 0. , ZSMAX , 0. , RS SMAX , 1, , ? 0 , XL AB EL , 10 , 
1 YLABSL.O) 

IEC =0 

DO 970 M=1,IPLTK 

RE AD (8 ) I NCOT(NI,N=l,L) 

IF (M.EO, IPLTK) I EC= 1 

970 CALL INFOPLT t IEC , L, X , 1. MOOT , 1 ,0. ,0. 1 0. »MDMAX, 1, , 1C , X2L , 1C , Y2L ,0) 
IEC =0 

IF (NNTP.EC.C) GO TO 1 
DO 990 M = ! , IPLTK 
I SYM=10 

DC 989 I = ! *NNTP 

READ (9) (ZZ(N) ,N=1,L) 

IF (M.EQ. IPLTK .AND. I.EO.NNTP) IEC =1 
I S YM= I SYN + l 

980 CALL INFOFLT ( I EC ,L ,X , l , ZZ , 1 ,0 . , 0 . , PTM I N , PTMAX , 1. , 10 ,X 2L ,20 , Y 3L , 

1 I SYM ) 

990 CONTINUE 
GO TO 1 
END 


33400000 
33500900 
33600000 
337C0CC0 
3 3 8000 00 
33900000 
34000000 
34100000 
34200000 
343000CC 
344CC0C0 
34500000 
346C0CC0 
3470CCCO 
348C00C0 
34900000 
35000000 
35 lOOOCO 
35200000 
35,300000 
354COCOC 
3550000.0 
35600000 
3570000C 
358 OOQCO 
359C00C0 
36000000 
36100000 
36230000 
363C00C0 
36400000 
36500'' CO 


Subroutine COLUMN .- Subroutine COLUMN calls the appropriate routines to com- 
pute the coefficient for the matrix solution and to solve the tridiagonal matrix for each 
column of temperatures. The flow chart for subroutine COLUMN is as follows: 




'CALL COLXO 
Compute 
coefficients 
i for column 1/ 


/CALL SOLMAT ’ 
/ Solve for 
\ TTF for 
\ column 1 
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The program listing for subroutine COLUMN is as follows: 


SUBROUTINE COLUMN 

SOLVES THE MATRIX COLUMN BY COLUMN FOR ONE ITERATION 

SOLVES M (NO. OF ROhSI SETS OF SIMULTANEOUS EQUATIONS N (NO. OF COLUMNS) 
TIMES THEN RETURNS TO MAIN PROGRAM TO TEST FOR CONVERGENCE 

COMMON /PICK/ A(lC,20t,AA(20),A8(10,20),ALPHA(20),B(20), 

2 BS1( 10,20 ,BS1B< 1C ,20 I, C( 10,20 ) ,CB( 10,20 ) ,CK( 10), CKETA( 10,20 ) , 

4 CKXI(10,20),C0ST(20),CP(10,20),0( 10 , 20) , DC ( 20 ) , 

6 OELESQ.DELETA, DELTA (20) ,D£LX I ,DELX ISQ , E ( 10 ,20 ) , E 1GHT3 , 

8 ELAM( 20) ,ETA( 10 ) , EXPG ,F ( 10 , 20 i ,GG,G I MACH ,H1 ( 10 ,20 ) ,H2( 10,20) , 

A H3(10,20),HC(20),HCOMB(2C),HE,HW(20),IFIRST,IROCOL,ITC,ITR,ITT, 

C ITT0.LM1 ,LM2 , MCDOT( 20 ) , MOOT ( 22) ,MSDOT(20) , P ID2, PRELOC ( 20) ,QC(20) , 

E QC1 ,QCNET(20) ,QCOMB(20) ,QR ( 20 ) , QRi, QS ( 20 ) , RNS ,RODPC ,ROPCPP, 

G RSS ( 22 ), RST02, S IG, S IGDP, SIGMA, SIGP, S INT( 2Q ) , SMI, SM2,TAU,TB, 

1 TT ( 10 ,20 ) ,TTF (10,20) , TWDELX I , TWOGI ,V(20) ,VB(10) ,X(22) ,XDXISQ, 

K XCDXI,Y(10,20),Z(20),ZB(10) 

COMMON /INPUTS/ DUMMY (1) , AEXP , ALCT AB ( 10 ) , TTALC( 10 ) , MAL PHC , NAL PHC , 

2 ALPHAT(IC) ,TALPHA(1C) ,MALPHA , NALPHA, ALST AB ( 10 ) ,TTALS( 10 ) , MAL PHS , 

4 NAL PHS, ASEXP, BETA,BEXP, BSEXP, CE, CKETATB ( 50 ) .TTCKETA ( 1 0) , 

6 ETATAB(5),NCKETA,NETA,CKXITAB(50),TTCKXI(10),XITAB(5) ,NCKXI, 

8 NXI , CORDS Y, CP DP ,CPP,CPTAB(10),TTABCP(10),MCP,NCP»DELTA0(20) , 

A DELTAU,DELTMIN,DTMAX,ELAMTB(28) ,TTELAM(7I , PELAMI4 ) , NELAM, 

C NPELAM,ENOTAU,EPSONE,EPSONEP, EPS0NPP,ERR0RT ,GAMBAR , C-AMI NF , 


36600000 

36700000 

36800000 

36900000 

37C00000 

37100000 

37200000 

37300000 

37400000 

37500001 

37600000 

37700000 

37800000 

37900000 

38000000 

381000C0 

38200000 

383C0000 

38400000 

38500000 

38600000 

38800001 

38800002 

38900000 
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E HCOHBTB { 28 ) » TTHCCMB( 7 ) » PHCOMB ( 4} , NHCOMB , NPHCOMB,HCT AB ( 2 8 S , 

G TTABHCJ7 8 »PHC(4) t NHC,NPHC »HETA3( 10) , TTABHE ( 10 ) s MFE «NHE ,HWTAB ( 15 ! » 
I TTABHW115) , MH4» *NHWj I ADJUST » I PLOT . L, MACHNO , MAX ITT »MDMAX , 

J MDCTO ( 20 ) * 

K MW02»MV(STR,NTP(7) ,PLTIME(15 ), PRAT {20! » PRFREQs PSEXP , PSTAGTBI 10 ) s 
M TTPSTAG(10)»MPSTAG ! NPSTAG»PTMAX,PTMIN,QCTAB(10) ,TTABQC ( 10 ) , MQC» 

N NCC i QRAT (20 ) • 

0 QRRAT ( 20 ) ,QRTAB(10t .TTABQR ( 10 ) , HQR # NQR ,R ( 20 ) , R I EXP, RNS I ,R0, ROOP , 

Q R0P,RS(20) ,RSSMAX,S, STEBOL,T( 10,20 ) , T AUO , TBTAB( 10) ,TTA8TB( 10) , 

S MTB,NTB,TDPRIPE, THETA (20) , TPR IME , XO, XOROE R, ZS (2C ) , ZSM AX 
REAL M DOT C » MOOT » MCOOT t MSDOT, Ml* STR » MWQ2 »MACHNO t MDMAX 
INTEGER S i SMI » SM2 
C COMPUTE COLUMN 1 
N1 =2 
N2 = SM1 

CALL COLXC (N1,N2) 

CALL SOLMAT ( A < 1 ,1 1 ,B ,C< 1 , 1 ) . Z ( 1 ) f V ( 1 ) ,DC,TTF( 1 , 1 ) , S ) 

C COMPUTE COLUMN 2 THRU LM1 
DO 300 N=2,LM1 
CALL COLMN (N1,N2,NI 

CALL SOLMAT ( A ( 1 , N ) , B ,C ( 1 , N ) ,Z (N ) , V (N) ,DC , TTF( 1 ,N ) ,S ) 

300 CONTINUE 
C COMPUTE COLUMN L 

CALL COLXL ( N1 » N2 » 

CALL SOLMAT ( A ( 1 ,L) , 8 ,C( 1 » L > , Z (L ) , V ( L ! ,DC, TTF( 1 ,L ) , S ) 

600 RETURN 
END 


390G0000 
39100000 
39200001 
39300000 
39400C00 
39500000 
39600000 
39700000 
39800000 
399C0000 
40000000 
40100000 
40200000 
40300000 
40400000 
405000 CO 
40600000 
40700000 
40800000 
40900000 
41000000 
41100000 
41200000 
41300000 
41400000 
41500000 
41600000 


Subroutine ROW .- Subroutine ROW calls the appropriate routines to compute the 
coefficients for the matrix solution and to solve the tridiagonal matrix for each row of 
temperatures. The flow chart for subroutine ROW is as follows: 


Q ROW ^ 


< CALL COLXO, 
COLMN, COLXL 
Compute 
coefficients 
for row 1 


CALL SOLMAT \ 
Solve for 
TTF for 
row 1 


CALL COLXOM, 
COLMNMN , COLXLM 
Compute 
coefficients 
for row 
2 thru S-l 
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The program listing for subroutine ROW is as follows: 


SUBROUTINE ROW 

SOLVES THE MATRIX ROW BY ROW FOR ONE ITERATION 

SOLVES N (NO® CF COLUMNS) SETS OF SIMULTANEOUS EQS. M(NC.CF ROWSI TIMES 
THEN RETURNS TO MAIN PROGRAM TO CHECK FOR CONVERGENCE 

COMMON /PICK/ A(1C,20),AA(20),AB(10,20),ALPHA(20),B(20), 

2 BS1(10,2C ),BS1B(1C,2C ),C( 10,20) ,C8UC,20) ,CK<10) ,CKETA( 1C ,20 ) , 

4 CKXI(iO,20) .COST (20) ,CP (1C , 20 ) ,D ( 10 ,20 ) , DC( 20 ) , 

6 DELES Q,DFL ETA, DEL TA< 20) , DEL XI , OE L X I SQ ,E ( 1 0 ,20 ) , E IGHT3 , 

8 ELAM( 20) , ETA (10 ) , EXPG , F ( 10 , 20 ), GG.GIMACH.Hl (10,20 , HZ (10, 20) , 

A H3( 10,20) ,HC( 20) ,HCOMB< 20 ) , HE ,HW ( 20 1 , I F IRST , I POCCL , ITC, ITR , ITT, 

C ITTC.LM1 , LM2 , MC DOT ( 20 ) , MOOT ( 22 ) , MS DOT ( 20) ,PID2iPPELCC( 20 ,QC( 20) , 

E QC1 ,QCNE T ( 20 ) ,QCOM8(20) ,CR ( 20 ) , QR1 , QS < 20 ) , RNS , RODPC .ROPCPP, 

G RSS(22) ,RSTC2 ,SIG,S I GOP, SIGMA , SI GP, SINT ( 20 ) ,SM1 ,SM2 , TAU , TB , 

1 TT( 10 ,20) ,TTF ( 10,20) , TWDE LX I , TWOGI ,V(20 ) , VB ( 1C) , X ( 22 ) .XDXISQ, 

K XODX I , Y ( 10 , 20 ) , Z ( 20 ) , ZB ( 1 0 ) 

COMMON /INPUTS/ DUMMY! 1) , AEXP , ALC T AB ( 10 ) ,TTALC(10) , MAL PHC, NAL PHC , 

2 ALPHAT (1C ) » T ALPH A ( 1 0 ) ,M AL PH A, NAL PHA , AL STAB ( 10 ) , TTAL S( 1C) .MALPHS, 

4 NALPHS, ASEXP,BETA,'BEXP,BSEXP,CE,CKETATB(50 ) , TICK ETA (1 0 I , 

6 ETATAB(5),NCKETA,NETA,CKXITAB(50 ) ,TTCKXI( 10) ,XITAB( 5) ,NCKXI , 

8 NX I ,CORDSY,CPDP,CPP,CPTAB(10) , TT ABCP( 10 ) , MCP , NCP , CELT AO (20 ) , 

A DELT4(J,DfcLTMIN,DTMAX,ELAMTB(28l » TTEL AM( 7 ) , PEL AMI 4 ) , NE LAM , 

C NPELAM, ENDTAU.EPS ON E,tPSONEP,EPSONPP,ERRORT,GAMBAR,GAMINF, 

E HCOMB TB (28) , T THCOMB ( 7 ) , PHCOMB ( 4 ) , NHCCMB , NPHCCMB ,HCTAB ( 2 8 ) , 

G TTABHC17) , PHC ( 4 ) , NHC , NPHC , HET AB ( 10 I , TT ABHE ( 10 > , MFF , NHE , HWTAB ( 15), 

I TTABHWU? ) , MHW , NHW, I ADJUST, I P LOT , L , MACHNO , M AX I TT , MD M A X , 

J MDOTO ( 20 ) , 

K MW02, MWSTR,NTP( 7) ,PLTIME( 15) , PRAT ( 20 ) , PRFR EQ, PSEXP, PSTAGTBt 10 I , 

M TTPSTAG(IO) , MPSTAG,NPSTAG,PTMAX, PTM IN , QCT AB ( 1C ) , TTA8QC ( 1C ) , MQC , 

N NQC , QRAT ( 2C ) , 

0 QRR AT (20 ) , QRT AB ( 1 0 ) , TTABOR ( 1C ) , MQR ,NQR , R ( 20 ) , R I E XP , RNSI ,PO,RGOP, 

Q RCP,RS(2C ) ,RSSMAX,S,STEBOL,T( 10 ,20) , T AUO , T BT AB ( 1C J , TT ABTB ( I 0 ) , 

5 MTB,NTB,TDPRIME»TFETA(20) , T PR I ME , XO , XCRDE R , l S ( 2C ) rZSMAX 
REAL MDOTC, MOOT, MCDOT.MSDOT, MW STR.MW02, MACHNO, MDM AX 


417GOCCO 
4180CCCC 
4190CCC0 
42000000 
421 OOCOO 
422C00C0 
423CCCC0 
42430CCG 
425C0C CO 
4260CC01 
427C0CC0 
428COOOC 
42 90 COCO 
43COOCCO 
43 iC DC CO 
432COCOO 
43300000 
43400CC0 
43 5 COO CO 
436C0C.CC 
4370 00 00 
4' 3 9000 01 
43900CC2 
4400 OCCO 
441CCCC0 
442C00C0 
443000C1 
444C00C0 
44S0C0C0 
446C0C00 
4470CCC0 
44300GCO 
4490 OC CO 
45C00CC0 
451C GCOC 
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C COMPUTE R )W 1 

01 PENS ICN AN S ( 20 ) 9 AT tMP ( 20 ) * CTEMP(2C( 

INTEGER S M I ,S 
N1 =2 
N2 = L M l 

CALL COLXC (NT ,N2> 

00 3)0 N = 2 t L M I 

CALL CCLMN (N1,N2,N> 

3.10 CONTINUE 

CALL C0LXL(N1,N2) 

DU 320 N =1 tL 
ATEMP(N) = AB(1,N) 

321 CTEMP(N) = C B ( 1 » N ) 

CALL SOLMAT < ATSMP , B ,C TEMP , Z B (1) , VB ( 1 ) ,0C , ANS (II , L I 
DC 400 N=3 ,L 

4) 0 TTF ( 1 * N1 = ANSI N ) 

C COMPUTE ROW 2 THRU SMI 
DO fel l M =2 t SMI 
N1 =M 
N2 =M 

CALL COLXCM (,N1,N2) 

DO 500 N = 2 9 L M ' 

CALL COLMNPN ( N 1 , N2 f N) 

5) 0 CCNTINUE 

CALL CQLXLM <N*_ ,N2 ) 

DC 51) N = 1 9 L 
ATEMPtN) = AB(N,N) 

51) CTEMPINI = CB(M,N) 

CALL SOLMAT ( ATEMP , B , CTE MP , Z B ( M ) , V 8 ( M 1 , DC , ANS ( 1 ) , L I 
DO 590 N= 1 9 L 
590 TTF ( M t N ) = ANS ( N ) 

S)C CONTINUE 
C COMPUTE ROW S 

CALL CQLXC1 (Nj ,N2) 

DO 8)0 N=2 9 LMi 
CALL CCLMM (N3 »N2,N) 

3)0 CCNTINUE 

CALL COLXLlINi ,N2) 

DC 810 N=l f L 
ATEMP(N) = AB( S 9 N ) 

910 CT-EMPI N ) = CB(S,N) 

CALL SCLMAT ( ATEM P , B , CT EMP , ZB ( S ) , VB ( S t , DC , ANS (1> , L ) 
DO 89) N=l f L 
99C TTF(S 9 N)=ANS(N) 

9)0 RETURN 
END 


45200000 
453C0C00 
454C0CCC 
4550CCCC 
456CC0C) 
45700000 
458C0CC0 
45900000 
46CC0CC0 
461 COO 00 
462C0CC0 
463C00C0 
46400000 
4650 CO 00 
46600000 
467CCCOO 
46800000 
469CC000 
4 7 C 0 0 0 00 
471C0C CO 
47200000 
473C00C0 
47400 COO 
475CCC0C 
476CD00C 
477C0CC.0 
478C0C00 
47900000 
4801CCC0 
48100000 
482C0CC0 
483000C0 
484COCOO 
485CCC00 
486C0C00 
48700000 
488C0000 
48900000 
49CC0CT0 
491COOOO 
492C00C0 
493CC00C 
494C0000 
495CCCC 0 
496COCOO 
497C0000 


Subroutine COLXO .- Subroutine COLXO computes the coefficients of the tridiagonal 
matrix where £ = 0 and 0 = ?y = 1. The flow chart for subroutine COLXO is as 
follows: 
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The program listing for subroutine COLXO is as follows: 


SUBROUT I Ml COL XC < M , N2 ) 
compute; coef. for xi=c, column implicit 

IRUCUL = 1 COLUMN) IMPLICIT 

IROCOL = 2 RCW IMPLICIT 

COMMON /PICK/ A( 10 ,20) ,AA( 20) , A3< 10 ,20) » ALPHA<20 ) , R<20 ) , 

2 BS1 (10.2C ) ,BS 16 HC, 2C ),C( 10,20) ,CB (10,20 ,CK( IC) ,CKETA(1C ,20) , 

4 CKXKiJ, 2T), C OS T( 20), CP (1J, 20), 0(10, 20), OC (20), 

6 DEL ESQ* DFLET A , DEL TA ( 20) , DEL X I ,0ELXISQ,E(10,2C) , EIGHTS , 

8 LLAMI 23) , ETA ( 1 0 ) , EXPG.F (10 , 20 ) , GG.GIM ACH ,H1 ( 10 , 2C ) , H2U C, 20 I , 

A H3( 13,20) ,HC(2C),HCOMB(2C ),Hc,HW(20) , I F I RST , I POCC L , I TC , I TR , I TT, 

C ITTC.LM1 ,LM2, MCC0TI2C ) , MOOT (22) , MSDOT ( 20 ) , P 1 02 , PRFLOC l 20 ) ,QC,< 20 ) , 
E QCl ,QCNET(20) ,QCCMB<23) ,0R(20) ,QR1,QS(2C) ,RNS,RODPC ,ROPCPP, 

G RSS (2 2) , RST02 ,SIG,SIGQP,S IGMA , S I GP , S INT ( 23 ) , SMI , SM2 ,T AU , T8 , 

1 TT( 10,20) ,TTF (10,20) , TWDE LX I , TWOGI , V ( 20 ) , VB ( iC ) , X 1 22 ) , X CX ISO , 

K X0DXI,Y(iC,23),Z(20),ZB(10) 

COMMON /INPUTS/ DUMMY (11 , AE XP , ALC T AB ( 10 ) ,TT ALC U C ) , M AL PHC , NAL PHC , 

2 ALPHATdO.T ALPHA (10) ,M ALPHA , NAL PHA , AL STAB ( 10 ) , TTALSI 1C) , MAL PHS , 

4 NAL PHS, ASEXP, 8ETA,BEXP,BSEXP,CE,CKETATB(50 ) , TTCKETA ( i 0 ) , 

6 E TATAB t 5 ) ,NCK ETA , NETA ,CKX I TAB ( 50 ) ,TTCKXI(10) ,XITAB<5) ,NCKXI, 

8 NXI ,C3RDSY,CPCP,CPP,CPTABU0 ) ,TTABCPI 10 ! , MCP , NCP , DELT AO ( 23 ) , 

A DELTAU,QELTMIN,DTMAX,ELAMTB(28),TTELAM(7) , PEL AM ( 4 ) , NELAM , 

C NPELAM, ENDTAU.EPSONE, EPSCNEP, EPSONPP , ERRORT , GAMBAR , GAMI NF , 

E HCOMBTB ( 28 ) , TTHCCMB < 7 ) , PHC OMB ( 4) , NHCC MB , NPFCQMB ,P CT AB ( 2 8 > , 

G TTABHC (7 ) , PHC 1 4 ) , NHC, NPHC, HETA8( 10 ) ,TTABHEI 10 I , MHE ,NHE ,HWTAB( 15), 
I TTABHtJl 15 ) ,MHW,NHW, I ADJUST, I PLOT , L , MACHNO , MAX I TT , MOM AX , 

J MD0T0I20), 

K MW02, M«STR,NTP(7) , PLT IME (15 1, PRAT (20 I , PRFREQ, PS EXP, PSTAGTB< 1C) , 

M TTPSTAGUC ) ,MPSTAG,NPSTAG,PTMAX,PTMIN,QCTAB( 1C) ,TTAB6C(1C) ,MQC, 

N NCC , QRAT ( 20 I , 

0 QRRAT(20) .CRT AB ( 10) ,T TABOR! 10 ) , MQR , NQR ,R ( 20 ) ,RIEXP,RNSI ,RO,RODP, 

Q RCP , RS (2C ) ,RSSMAX,S, STEBCL.TI 13,20) ,TAUC,TBTAB( 10 ) , TT ABTB ( 1 0 ) , 

5 MTB,NTB,TDPRIME,THETA(20) , T PR IME ,XO,X ORDER, ZS (2C- ) , ZSMAX 
REAL MOOT C » MDOT , MCDOT , MS CCT, Min STR »MW02 , MACHNO , MDMA X 
INTEGER S , SMI , SM2 


498CCC0Q 
499C OCC-U 
500C0C00 
507 COCCO 
502CCf CO 
E02C000C 
5G4CCCC0 
505C3C0O 
5U6COOOO 
5C7COOC1 
5C8CCCC0 
so sc croc 

5 ICC CC 00 
5 1 1 C C 0 00 
512COOOO 
5 1 30CC00 
5J 400000 
51 5C0C00 
51600COC 
51 7C COCO 
516C0C00 
52CC00C1 
52C000C2 
5210000C 
522C0CC0 
523C000G 
5240000) 
52500000 
52 EOCCOO 

52 70 OC CO 
528CCOOO 
52900000 

53 30G00C 
53100CC0 
53200000 
5 3 30 CC 00 
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STATION ( 1,11 XI=C , ET A =’3 
DO 60 1 = 1 , SMI 

60 CK ( I ) = (CKETAt 1 , 11 + CKETA ( 1 + 1 , 11 1 / 2.0 
DELDE = DELTA!! I * DELETA 
PART 2 = HI ( 1 , 1 ) **2 * XDXISO 

PAPTl=RQDPC 
H 1 R = Hid ,11 * R ( 1 ) 

FF=CKXI( 1 , 1 )*( 2 . 0 -C 0 RDSY)/( 2 . 0 *PART 2 ) 

G=RO*CP( 1 , 1 ) / 0 ELTAU- 2 . 0 *PART 1 /H 1 R+ 8 . 0 *PARTI/ ( 3 . 0 *DELCE ) 

H = 1 , 3 /( H 2 ( 1 , i ) **2 * DELTA ( 1 ) ** 2 ) 

SC= H / ( 3 » 0 * DELESO) 

FPT 4 =S IGDP* ( 2 » 0 /(HlR*H 2 (l , 1 ) ** 2 ) - E IGH T 3 /DELDE ) 

EPTB= EPT 4 *TB 

EPT 4 = £P T a *T( 1 , 1) **3 

BSAVE = G 

GO TO ( 70 » 80 ), IROCC'L 
70 CONTINUE 

A ( 1 , 1 ) = 0.0 

BSK 1 , 1 ) = -SC* 9 .C *CK ( 1 ) 

C ( 1 , 1 ) = SC 4 ( 9 .G *CK< 1 ) + CK ( 2 ) ) 

Z(ll = -SC * CK ( 2 ) 

B( 1 ) = 8 SI ( 1 , 1 ) - ESAVE + EPT 4 
IF ( IF IRST.EO.C > GO TO 80 

78 DC ( 1 ) = ( - RSAVE-BS 1 8( 1 * 1 ) ) *T ( 1 , 1 ) - C B ( 1 , 1 ) *T ( 1 , 2 )- ZE( 1 )*T( 1 , 3 ) 

1 + EPTB 
GO TO =) 

30 FP=FF 

BS 1 B( 1 , 1 )= - 7 . 0 * FP 
CB ( 1 , 1 ) = 6.0 *FP 
ZB ( 1 ) = -FP 

IF ( 1 FIRST.EQ . 0 I GO TO 78 
36 8 ( 1 ) = BS 1 B( 1 , 1 )— BSAVE + EPT 4 

DC ( 1 I =( — BSAVE -BSK! , 1 ) )*T(l,i) -C ( 1 , 1 ) *T ( 2 , 1 ) - Z( 1 )*T( 3 , 1 > 

1 + EPTB 

93 GO TO ( 101 , 600 ) , IROCOL 

ST AT I ON ( M, 1 ) , xi = 0 , ETA LESS THAN 1 , GREATER THAN 0 

ENTRY COL X CM 
131 DO 20 0 M=N 1 ,N 2 

DELDE = OcLT A ( 1 ) * DEL ETA 
MP 1 =M + 1 
MM 1 =M -1 

P 817 = 8 . 0 * DELTA! 2 ) - DELTAI 3 ) - 7 . 0 * DELTA(l) 

PAPT 2 = H! ( M , i ) **2 * XDXISO 

CORO= ( 2 .C-CORDSY 1 / 2 .C 
FF= CKXI(M, 1 )*C 0 RD/PART 2 
G = RO *CF(M, 1 ) /DELTAU 
SC = 1.0 / ( H 2 ( M , 1 ) * DELOE ** 2 > 

H = FF* P 817 / ( 2 . 0 * DELDE) *ETA(M) 

P = CKETA(M, 1 )/(H 2 (M, 11**2 *Hl(M,i>* R(i) * DELDE) 

BSAVE =G 

GO TO ( 17 C , 180 ) , IROCOL 
170 CONTINUE 

U= ETA(M) *MD0T( 1 ) * CP ( M , 1 ) / ( 2. 0*DELT A ( 1 ) * DELETA) 

A ( M , 1 ) = F -P + SC* CK(MMl) +U 
BS 1 ( M, 1 ) = SC * (-CK(MMl) - CK ( M ) ) 

C(M,l)= -F + P + SC* CK ( M ) -U 
B ( M I = 8 SH M ,1 ) - BSAVE 
IF (IFIRST.EQvO ) GO TO 180 

173 OC ( M ) = (-BSAVE -BS1B(M,1 ) )* T ( M, 1 1 -ZB ( M ) *T( M, 2 )-CB( M , 1 ) *T ( M ,2 ) 

GO TO 200 
130 ZB ( M ) = -FF 

CB ( M , 1 ) = 8.0 * FF 
BS 1 B(M, 1 )= - 7 . 0 *FF 
IF l IF IRST.EQ.c I GO TO 178 
190 B ( 1 ) = BSI B(M , 1 ) - BSAVE 

DC ( 1 I = (-BSAVE - BS 1 (M, 1 ) )*T(M, 1 )-A(M, 1 )*T(MM 1 , 1 )-C(M, 1 )*T(MP 1 , 3 ) 
230 CONTINUE 

GC TO ( 202 , 603 ) , IROCOL 


534GCCC0 
535CCOCC 
536COOOC 
537C0C OC 
538CCCCO 
53 9C 0000 
54CC OCOO 
541 CCOCC 
542C00C0 
543C.OOOO 
54hCC0CC 
54500000 
5460 000-0 
5470DC00 
548CCCC0 
549000CC 
550C00C0 
55100000 
552CCCOO 
55200000 
5540C00C 
555COOOO 
556CC00C 
557CC0C0 
5580C0C0 
559C OOCO 
5600 0000 
561CCC00 
56 20 CO OC 
563CG0C0 
56400900 
5 6 5 0 0 0 CC 
5660 0000 
56 7C0000 
56800000 
5 6 9C OOCO 
57000000 
57100000 
57200000 
57300000 
574CC000 
575C0000 
57690000 
577CCC00 
57800000 
57900000 
58000000 
5813CC00 
58200000 
583CCOCO 
58400000 
58500000 
58600000 
587CCCOO 
58800000 
589000C0 
59CC0000 
591C0GCC 
5920000C 
59300000 
59400000 
595C0000 
59600000 
59700000 
598GQCOO 
5990C0CC 
60000000 
60100000 
602C000G 
60300300 
60400000 
60500000 
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STATION (S,1S ,XI=0 , ETA=1 

ENTRY COL XC1 

202 CORD=( 2,0~CCRDSY)/2a0 

FF=CKXI(S,l)*CORD/Hl( S,l)**2 
DEL0E=DELTA(1)*DELETA 
P =FF/XDXISQ 

H = 1. 0/ ( H2 ( S , 1 ) **2 *3,0* DEL0E**2) 

G = RO* CPIS, 1) / DELTAU 
SC = -9.0 * CK (SMI I * H 
BSAVE = G 

GO TO 1270,280) , IROCOL 
270 CONTINUE 

XX=CP( S,1 )*MDOT( 1 )/(2.0#0ELTA(l)*DELETA) 

VI 1 > = -CK(SM2)*H - XX 

A( S , 1 1 = -SC + CKI SM2 )*H + 4.0*XX 

DR=P*P817/CKETA(S,1) *H2(S,1> 

DD = DR - 2.0/(Hl(S,l)*R(l )*H2(S,1))-EIGHT3/ 

1 ( H2 ( S, 1 ) *DELDE) 

DDQS=DD*QS(1) 

BSU S,1 )=DD*SIG*T( S,1 ) **3 +SC-3.0*XX 
BIS) = 8SI I S , 1 ) -BSAVE 
IF (IFIRST.EQ.O ) GO TO 280 

278 DC I S I = DCQS + (-BSAVE - BS1BI S , 1 ) ) *T ( S , 1 »- CB I S, 1 1 *T IS, 2 ) 

l - ZBIS) *T I S , 3 ) 

GO TO 600 
230 CBIS.i ) = 8.C*P 
ZBIS) = -F 
BS1B I S , 1 ) = -7.0*P 
IF (IFIRST.EQ.O I GO TO 278 
290 BIl) = 8S1 B I S , 1 1 - BSAVE 

DC 1 1 ) = (-BSAVE - BS1(S,1 ) )+T(S,l) - V ( 1 )* 

1 T( SM2, 1 ) - A(S,1) *T ( SMI , 1 ) +DDQS 
600 RETURN 

2n FORMAT (7E16.7) 

END 


606C00CC 
6070C0C0 
608CC000 
60900000 
610C0000 
61100000 
612C0000 
61300000 
61A00000 
61500000 
61 6C0000 
61700000 
61 800000 
61900000 
62000000 
621C0C00 
622COCOO 
62 3C 0000 
62400000 
625C00C0 
62600000 
627C0000 
62800000 
62900000 
63000000 
631C0000 
63200000 
63300000 
63400000 
63500000 
636C0C00 
637G0C 00 
63800000 
63900000 
64COOCOO 
641COOCO 
64200000 


Subroutine COLMN .- Subroutine COLMN computes the coefficients of the tridiagonal 
matrix where 0 < £ < 1 and 0 ^ r] i 1. The flow chart for subroutine COLMN is as 
follows: 
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The program listing for subroutine COLMN is as follows: 


SUBROUTINE CCLMN (NI,N?,N) 

I ROCOL = 1 COLUMN IMPLICIT 

IROCOL = 2 ROW IMPLICIT 

COMMCN /PICK/ AI 10,20) »AA( 20) , A, 3{ 10,20 »ALPHA(2C) ,B(29 > , 

2 BS1 1 10 , 20 ) ,BS1P. (10,20 ) ,C( 10,20) »CB( 10,20 ,CK( 10) ,CKETA( 1C, 20 I , 

A CKXI( 10,20) ,CCST( 20 ) ,CP< 10,20) ,0(10,20) ,DC( 20) , 

6 0ELESQ,0ELETA,DELTA(20) ,DE LX I , OELX I SC , E { 10 , 20 ) , E IGHT3 , 

8 EL AMI 20) ,ETA( 10) , EXPG ,F ( 1C , 20 ) ,GG,GIMACH ,H1 ( 1C , 20 ) ,H2 (10,20) , 

A H3 ( 10 , 20 ) , HC t 20 ) , HCOMBI 20 ) , t-E , HW ( 20 ) , I P IRST , I P0C0.L , I TC , ITR , I T T , 

C ITTG,LMl,LM2,MC0CT(2O .MOOT (22 1 ,MS00T(2C( , PIOZ , PPELOC ( 20 ) , OC ( 20 ) , 
E QC1 , QCNET ( 20 ) , QCCMB ( 20 ) ,QR ( 20 ) , QR1 , QS ( 20 I ,RNS ,RCDPC .ROPCPP, 

G R SSI 2 2) ,RST02,SIG,SIGDP, SIGMA, SIGP,SINT(20), SMI, SM2,TAU,TB, 

1 TTI 10 ,20 ) , TTP 1 10 ,20 ) , TWDE LX I ,TWOGI,V(20),VB(1C) ,X(2 2) .XOXISQ, 

K XODXl ,Y( 10,23) ,2(20) ,ZB(10) 

COMMON /INPUTS/ DUMMY! 1 ) ,AEXP,ALCTAB(1Q) ,TTALC(10) ,MAL PHC,NALPHC , 

2 ALPHAT CIO) ,T ALPHA (1C I , MAL PH A , NAL PHA , AL ST AB ( 10 ) , TT AL S I 1C ),MALPHS, 

4 NALPHS, ASEXP, BETA , BE XP , BSE X P , CE , CKE TATB ( 5 0 ) .TTCKETAIIO) , 

6 ETATABI5 ) ,NCK ET A , NET A ,CKX I T AB ( 50 ) ,TTCKXI ( 10 ) ,XITAB( 51 ,NCKXI , 

8 NXI ,C0R0SY,CPCP,CPP,CPTAB<10> ,TTABCP(10) , MC P , NCP , DE LT AC ( 20 > , 

A 0ELTAIJ,DELTMIN,DTMAX,ELAMTB(28) .TTELAMI7) ,PELAM(4) » NE L AM , 

C NPELAM.ENDTAU ,EPSCNE , EPSONE P , EPSCNPP »ERRORT, GAME/ R , GAMINE, 

E HC0M8TB(28) ,TTHCCMB(7 ), PHC0MB(4I , NHCOMB .NPHCOMB ,HCTAB ( 2 8) , 

G TTABHC(7) ,PHC (41 , NHC , NP HC , HET AB ( 10 ) , TT ABHE 1 1 0 ) , ME E , NHE , HWT AB ( 15 ) , 
I TTABHW( 15) ,MHW,NHW, I AD JUST , I PLOT , L , M ACHNO , MA X I TT , MDMA X , 

J MDGTD ( 2 3 ) , 

K MW02, MWSTfi,NTP(7) .PLTIME ( 15 ) , PR AT (20) , PRFREC, PSEXP , PSTAGTBI 1 0 ) , 

M TTPSTAGt 10 ) , M PSTAG, NPST AG, PTM AX , PTM IN , OCT AB 1 10 ) , T TABOC ( 10 ) ,MQC, 

N NOC , QRAT I 20 ) , 

C CRRAK20 ) , CRT ABO 01 ,TTABQR( 10 I , MQR , NQR , R ( 20 1 , R I EXP , RN SI ,RO,RODP, 
0 R0P,RS(2C),RSSMAX,S,STEB0L,T( 10,20) , T AUC , TBTAB 1 10 ) , TT ABTB I 10 ) , 

5 MTB,NTB,T0PRIME,THETA(20),TPRIME,X0,X0RDER,ZS(2C),ZSMAX 
REAL MD0TC,MD0T,MCD0T,MSD0T,MWSTR,MW02,MACHNC,MDMAX 
DIMENSION DDQS ( 20 ) , ODQ SR ( 20 ) 

INTEGER S , SMI , SM2 

STATION ( 1 , N) XI GREATER THAN 0,LESS THAN 1 ET A=0 

231 FORMAT (7E18.7) 

NM 1 = N-l 
NP1 = N+l 

E32N=(H1( 2,N)+HT(1,N) >*{H2 (2 ,N)+H3(1 ,N) ) * ( CKETA I 2 ,M +CKET A ( 1 ,N) )/ 

1 !4.*(H2(2 ,N)+H?( 1 ,N) ) I 

E52N = ( H1(3,N) + H1 (2,N) )*(H3(3 ,NKH3(2,N) ) * < CKETA( 3 ,N ) +CKETAI 2 , N ) 1 / 
l (4«*(H2(3,N)+H2(2,N) ) I 
VV=1 • 3 / ( 3 » C * DELTA! N) **2 * DELESQ ) 

P1NP1=(H3I1,NP1 )+H3( 1,N) ) / I HI ( 1,NP1)+H1( 1,NI ) 

PI NM 1= (H3 1 1 , NM 1 ) +H31 1 , N) ) / ( HI ( 1 , NM1 ) +H1 1 1 ,N ) ) 


64 3 COO 00 
64400000 
64500000 
64600000 
647G0CC0 
64800000 
64930000 
65QOOD n O 
651C0C01 
65200000 
65300000 
6540G00C 
655C0000 
65630000 
65700000 
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65900C00 
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66 1C 0000 
662CC000 
66400C01 
66400002 
66500000 
6660000C 
667C0000 
668C0CC1 
669C0C00 
67000000 
671G0CC0 
67200000 
67300000 
674000C0 
675C00C0 
6 76 OC 000 
677C0000 
67800000 
679COCOO 
68DOOOCO 
681 00000 
6 82000C0 
68200CCO 
68400000 
68500000 
686C0Q0C 
68700000 
68800000 
68900000 
690C0CT0 
69 1 COC 00 


32 






o o n o 


W= Hl(l,M*H5U,M* CELETA *OELTA(NI *8»C 
GIN = Hid.N)* H2(1,M * H ? { 1 , N ) * RO *CP(1,N) 
YY=(-VV*W*R0DPC-G1N/DELTAU) 

C PT4 = -VV *W * S I GDP 
EPTR= EPT A * Tfi 
EPT4 = EPT4 * T(1,N)**3 
BSAVE = YY 

GC TO ( 170(18'' ) > IFCCOL 
170 CONTINUE 

BSl(lfN) * “VV* 9.0* F 32N 

C ( 1 , N I = VV *(9.0* E22N + F52M 

Z ( N ) = -VV * E 5 2N 

8(1)= 8S1 ( 1 *N) + PSAVE + EPT4 

IF ( IF IRST.EO. C ) GO TO 180 

178 DC ( 1 ) = (BSAVE - BSi 8 ( 1 , N )) *T ( 1 , N ) - A8 ( 1 , N ) *T ( 1 , KM 1 ) -CB ( 1 , N ) * 

1 Ttl.NPl) + EP Tfi 
GO TO 200 

180 D1NP1=(H2(1 ,NP1 l+H2(l ,NI )*(H3( i ,NP1 ) +H3 ( 1 , N ) ) * ( CKX I ( 1 * NP1 l+CKX 1(1 , 
IN) )/ (4, *X CXI SO* (FI (1,NPU + H1( 1 ,N) ) ) 

D1NM1 = ( H2 <1 » KM1 ) *H2( 1 » N) )*(H3(1«NM1)+H3(1»N) )*(CKXI ( 1 , NM1 ) +CKX I ( 1 , 
IN) )/<4,*XCXIS'3*(Hl(i,NMil + Hl(l,NI I I 
ABU. ,N) =D1KM1 
BSiB(l,'O=-01NPl- D1NM1 
C8( 1 , N I -0 • KPl 

IF ( IF IRST.EQ. 0 ) GO TO 178 
190 R(M = BSiPd.N) + BSAVE + EPT4 

DC(N) = ( PSAVF -BS1(1,N) >*T(l,N) -C ( 1 ,N ) *T ( 2 ,N ) - Z(K)*T(2,N) 

1 + EPT3 
2)0 CONTINUE 

GC TO (202,800 1 , I ROC GE- 
STATION ( M , N ) XT GREATER THAN C, LESS THAN 1 
ETA GREATER THAN C, LESS THAN 1 

ENTRY COLMKMN 

NP1=N+1 

KMi=N-l 

2)2 DO 400 M=M ,N2 
MMi = M-l 
MPI = M+ 1 

VV= 1,0 / (DFLTA(NI**2 * DELE SO) 

XX = cTA(f')*AA(N)/(D£LTA(N)* DELE SC ) 

G = HI ( M* K) * H2(M,NI * H3(M,N) *RO * CP(M,N) * 

EMM12N = (H1 ( MMi ,N)+Hl <N,N) ) *(H3(MMi ,N) +H3(M,N) ) * ( C K ET A ( MMI , N ) + CKET A 
1 < M,N) )/(4, *<H2(NMi,N)+H2(M,NI))*VV 
EMP1 2N = (H1 (MPI ,K)+H1(M,NH*(H3(MP1,N)+H3(M,N) )*(CKETA(MP1,N)+CKETA 
1(M»N))/(4,*(H2(MP1,N)+H2(M,N)))*VV 
DKM12N = (H?( MMI , N I +H7 ( M , M ) ) * ( H3 ( MM1 , N I +H3 ( M , N ) ) * ( CK X I ( MM 1 , N ) +CK X I ( M 
1 M 1 , N ) ) / ( 4 „ * ( H 1 ( MM 1 * N ) + HI ( M , N ) ) ) 

DMP12N = (H2( MPi ,N)+H2( M,N) ) *(H3 ( MPI ,N)+H3(M,N) ) *< CKXI (MF1 ,N) + CKXI(M 
1 PI ,N))/(4i*(Hl(MPi,N)+Hl(M,N) ) I 
FNM12N=XX*DMM12N*AA(N)*(ETA(MM1I+ETA(M) ) / ( DELTA ( N ) *2. ) 

FMPT 2N=XX*DMP12N*AA(N)*(ETA( MP1)+ETA(M) ) / ( DELT A ( N ) *2, ) 

W = *. 0 * XC * CELXI * DEI ETA 

DENOM= 4. C * ( DELTA? NM1 ) + DELTA ( Nl) *( HUM, MU) + HI (M,N)I 
FMNM12= (H3(M»NM1)+H3(M,N) ) * ( H2 { M, NM1 ) +H2 ( M ,N ) ) * ( CKX I ( M, NM 1 ) 

1+ C K X I ( M , M | * ( A A { NM1) + A A ( N))* ETA ( M I /DEKOM 

DENOM= +.C* (DELTA) NP1 ) +DELT A( N > ) * < HI ( M , NP 3 > +HJ ( M, N > I 
FMNP 1 2 = (H3(M,NP1)+H3(M,N) )*(H2(M,NP1) +H2 ( M , N ) ) * ( CKX I ( M , NPl ) 

1 +CKXI ( M , M )*( AA( NP1I+AAI N ) ) *E T A ( M ) /DE NOM 
D1 = (FMNPI2*(T{MP1,NP1)-T(MM1,NPII+T(MP1,N)-T(MM1,N>)-FMNM12* 

1 (T( MPI ,N 1-TtMMl ,N )+T( MPI, NM1 )-T( MM1,NM1 1 ) ) ZW 
02 = ETA(M *AA(N»* CKXI (M,NI* (H2(MP1,N)* F3(MP1,M* (T(MP1,NP1) 

1 - T(MP1,KM3. ) ) /H1(MP1,N) - H2(MM1,N) * H3(MM1,N)* (T(MM1,NP1) 

2 - T(MM1,NM1) )/Hl(MMl,N) ) / (DELT A ( N I * W) 

DS =01 +02 - G *T(M,N )/ DELTAU 

RSAVE = G/DELTAU 

GO TO (37C.33C) , IPOCCL 
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694000C0 
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3 70 CONTINUE 

HMN = ETA ( R ) * MDOT(N)/(DELTA(N) * R0) 

YY = G * FRN/ ( 2.0* DELETA) 

A(M,NI = t RR12 N + FMM12N + YY 

BS 1 ( M, N ) = -EMM12N - EMP12N -FMPI2N - FRR12N 

C(M,N) = FRP12N + FMPI2N - YY 

B ( M » = BSIIM.N) - BSAVE 

IF (IFIRST.EQ.O ) GO TO 230 

378 DC ( M ) = DS- BS 1 8 ( M ,N ) * T (M,N) -AB I M, N) *T ( M ,NM1 ) -C B I M, N 1 *T ( M, N PI ) 
GC TO 430 

380 DMNMi 2= ( H2(M*NR1)+H2(M»N) ) * I H? ( M, NM1 1 +H3 ( M , N ) ) * ( CKX I t M * NM1 ) +CK X I ( M 
1,N)>/(4.*(H1(M,NM1MH1(M,N> J ) 

AB(M,N)=DMNM12/XDXISQ 

DMNP12=(H2(R,NPi ) +F2 t M , N I) * ( H3 { M, NP1 )+H 3 ( M , N ) ) * ( CK X It M , NP I ) + C K X I ( M 
l,N))/(4.*IHl(M,NPl)-t-Hl(M,N))) 

CB(RI,N)=DRNP12/XDXISQ 
BS1B(M,N)= -AB(M,N) - CB(M,N) 

IF (IFIRST.EQ.O I GO TO 378 
3°0 BIN) = 8 S3 B I M * N ) - BSAVE 

DC ( N ) = DS - BS1(R,N)*T(M,N) - A ( M, N ) *T ( MR 1 , N ) - C ( R , N ) *T ( RP1 , N ) 

430 CONTINUE 

GO TO (431,830) , IROCOL 

STATION (S,N) XI GREATER THAN C, LESS THAN 1 , ETA =1 

ENTRY C0LMN1 
NP1=N+ 1 
NMl=N-l 

431 F1H3 = HKS.NI* H3(S,N) 

XX= 3,3 * 0ELTA(N)**2 * DELESC 

U = AA ( N ) / (3.0 *DELESC* DELTA(N) ) 

G = hiH3 *H2(S,N) * RO *CP(S,N) 

PART = AA(N)/(D2LTA(M*4.0*0ELETA*DELXI*XD I 
SST= F3(S,M*CKXI(S,N)*3./H1(S,N) 

DS=P ART* ( SST*T ( S »NP1)-SST*T( S , NM1 ) 

1 -4.3*F2(SM3,N)*H3(SM1,N)*CKXI( SM 1 ,N > * ( T ( SMI , NP1 ) -T ( S M 1 ,NM1) )/ 

2HKSM1 , N) + H2 ( SM2 »N)*H3 ( SM2 , N ) *CKX I ( SM2 , N ) * ( T ( SM2 , NP1 ) -T (SM2 , N Ml I I 
3 /HI ( SM2 , N ) ) 

cSM32N = t H 1 ( SRI , M +H1 ( SM2 , N ) ) * ( H3 ( SMI , N) +H3 ( SM2 , N ) ) * ( CKET A( SMI , N ) + 

1 CK ET A ( SM2 , N ) 1/(4. *(H2 (SMI, N)+H2(SM2,N) )*XX) 

ESM12N = (H: (SMI ,N)+H1(S,N I ) *( H3 ( SMI , N ) + H3 ( S , N ) )*(CKFTA( SM1,N)+CKET4 
1(S,N) I / ( 4. * ( H2 ( SMI ,N)+H2(S,N) ) *XX ) *9. , 

DSM1 2N= ( H2 (SMI ,N)+H2(S,N))*( H3 ( SMI , N ) +H3 < S , N ) ) *( CKXI ( SM2 ,N ) +CK XI ( S 
1 , N ) ) / (4. C* ( HI ( SM 1 » N ) + HI ( S , N ) I ) 

FSM12NNDS R 12N* AA (N I* ( ETA ( SMI ) + ETA < S I ) / ( DEL TA ( N )*2. ) *9. *U 
DSM3ZN=(H2(SM? ,N I +H2 ( SMI , M ) * ( H3 ( SR!2 ,N) +H3 ( SMI , N ) ) * ( CK XI ( SM2 , N ) + 

1 CKX K SM1,N) )/ (4.*(H1(SM2,N) + H1( SM1,N) ) ) 

FSM32N=DSR32N*AA(N)*(ETA( SM2 )+ETA(SMl) ) / ( DELT A ( N ) *2. )*U 
BSAVE = G/DELTAU 
GO TO (570,530 .IROCOL 
570 CONTINUE 

YY=G*MOOT ! M / ( R0*2. 0*DELTA ( N ) *DELET A ) 

V ( N ) = -ESR32N - FSM32N -YY 

A ( S , N ) = E SMI 2 N + ESM32N + FSM12N + FSM32N + 4.0*YY 
0D=8.*H1H?*DELTA(N!*DELETA/XX »■ 8. *U* 

1H2(S,N)*DELTA(N)*F(S,N )*OELETA/CKETA(S,N) 

DDQS(N)=DD*QS(N) 

BS 1 ( S, NI = -OD*S IG*T ( S , N ) **2-E SM12N-FSM1 2N- 3. 0*Y Y 
BIS) = BS I ( S , N ) - BSAVE 
IF (IFIRST.EQ.O I GO TO E8C 

578 DC ( S ) =-DD CS ( N ) + DS+ (-BS AVE-BS1 B ( S , N ) I *T ( S , N )- AB ( S , N ) *T ( S , NM 1 ) 

1 -CB(S,N)*T(S»NP1)+ ODQSRINI 
GO TO 530 

5 33 DSNM12=(H2 (S,NR1 1+H2I S,N) ) *(H3( S , NM1 ) +H3 ( S , N > )*(CKXI ( S , NM1 1 +CK X I ( S 
l,N)J/(4,*(HI(S,NMi)+Hi(S,M> *XDX I SQ ) 

0SNP12 = (H2(S,NPU+H2(S,N))*(H3(S,NP1)+H3(S,N) )*(CKXI(S ,NP1 l+CKXUS 
1 ,N) ) /(4,*(H1 (S ,NP1)+H1 (S,N ) ) *XDXISQ) 

DEN0M=4« C* ( DEL TA ( NM1 ) *DE LT A ( N ) ) *( H2 ( S , N Ml I +H1 ( S , M ) 

F SNM12 = (F3(S,NM1)+H?(S,N) ) * ( H2 ( S , NM1 ) +H2 ( S , N ) )* (CKXKS.NMl) 
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1 +CKXI (S,M )*< A A ( NMD+AAI NH/DENCM 
DENCM=4.0M0ELTA( NP1H-DELTA1 N )) *( HI { S » NP1 ) +H1( S ,N M 
FSNP12= (H3(S,NPl)+H3(S,Nn*(H2(S,NPl)+H2(S,N) ) *< CKX l ( S, NP1 ) 
1 +CKXI (S»N) )*< AA( NPI)+AA< N > I / DENQM 
DEN0M=2»0*XC*DELX I 

CSN= OELTAIM *H2 (S»N )/ (CKETAI S»N) *DENCMi 

OSNPl = OELTA(NPX)* H2 ( S , NP1 ) / ( CKETA< S ,NP1 ) * DENCM) 

QSNM1= DELTA(NMX)* H2 < S , NM I ) / ( CKE TA ( S , NM1 ) * DENOM) 

DDQSR ( N > = FSNPX2* < QSNPX*QS ( NP1 > + QSN*QS<Nl I-FSNM12* 

1 (QSN*QS<N)+ QSNM 1*GS ( NM1 ) ) 

AB ( S ,N)=DSNM12-FSNM12*SIG*QSNMX*T<S,NM1 1**3 
CB(S,N >=DSNP12+FSNPX2*QSNPX*SIG*T( S,NP1!**3 

BSXB ( S » N ) =-DSNP12-DSNM12 +S IG*T ( S » N ) **3 *QSN * ( F SNPX 2-FSNMX2 ) 
IF ( IFIRST.EQbO ) GO TO 578 
590 B(N)=8SlB(StM-PSAVE 

DC ( N I = -CD QS ( N ) +(-BSAVE-BSi(S*N) ) *T ( S » N I + DS 
X-A(S.N)*T<SMX,N)-V<N)*T(SM2,N>+ DDQSR (N) 

600 CONTINUE 
800 RETURN 
' END 


82900000 
83CCOCCO 
83XOCOCO 
832C00C 0 
83300C00 
8 3400000 
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83 6C00C0 
837CC0C0 
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843 C 0000 
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Subroutine COLXL .- Subroutine COLXL computes the coefficients of the tridiagonal 
matrix where £ = 1 and 0 = t] = 1. The flow chart for subroutine COLXL is as follows: 



35 










non o o o o r> 


The program listing for subroutine COLXL is as follows: 


SUBROUT IMF C0^XL(N1,N2) 

COMPUTES COEF, FOR XI=1 ( X=L ) COLUMN IMPLICIT 
IROCGL = 1 COLUMN IMPLICIT 

IROCOL = 2 ROW IMPLICIT 

COMMON / PICK/ AC 10,201 ,AA( 20) , A8< 10 ,20 , AL PH A( 20 ) , B ( 20 ) , 

2 BS1Q3.2C ),BS18(10,20 I ,C< 10 , 20 1 , CBI 10, 20 » , CK ( 10 ) ,CKETA( 1C ,20 I , 

4 CKXK 10,20) .CCSTC20I , CP (10, 20) ,D (10,201 ,DC(20), 

6 CELESQ,0ELETA,DELTA(20) , DELX I , DELX I SQ , E (10,20) ,E1GHT3, 

8 EL AMI 20) ,ETA< 10 ) , EXPG ,F (10 , 20 I ,GG,GI MACH ,H1 ( 10 ,20 ) »H2 ( 10, 20 ) , 

A H3 1 10 , 20 ) ,HC( 20 ) , HCOMBC 2C ) »HE»HW( 20) , IF IRST » l ROCCL » ITC , ITR , I TT, 

C ITTC,LM1,LM2,MC0CT(20) , MOOT (22) ,MS00T(20) , P 102 , PP ELCC ( 20 ) , OC ( 2C ) , 
E QC1 ,QCNET (20 ) ,QCCMB(20) ,0R ( 2C ) , QR1 , QS( 20 ) ,PNS,RODPC ,ROPCPP, 

G RSSC2 2) ,RST02 , S I G, SI GOP ,S IGMA ,S IGP , S INT ( 20 ) ,SM1 ,SM2,TAU,TB, 

1 TT(10,20) ,TTF(1G,20),TWOELXI , TWOGI , V( 20 ) , VB ( 1C) , X (22 ) .XOXISQ, 

K XCDXI ,Y( 10,20 ) ,Z( 20) , ZB(10) 

COMMON /INPUTS/ DUMMY ( 1 ) , At XP , ALC TAB ( 10 ) ,TTALC(10) , M AL PHC , NAL PHC , 

2 ALPHAT (1C) ,T ALPHA (10) , M AL PH A, NAL PHA , ALST AB ( 10 ) ,TTALS( 10 ),MALPHS, 

A NALPHS,ASEXP,BETA,BEXP,BSEXP,C£,CKETATB(50) .TTCKETAtlO) , 

6 E TAT ABC 5) , NCK E TA ,NETA , CKX IT AB ( 50 I .TTCKXI ( 10) ,XITAB( 5) ,NCKXI , 

8 N XI , CORD S Y ,C p DP,CPP,CPTAB(10) ,TTABCP(10) , MC P , NCP , OE LT AO ( 20 ) , 

A OELTAU,DELTMIN,DTMAX, ELAMTB(28), TTELAM(7) , PEL AM ( * ) , NE LAM, 

C NPELAM,ENDTAU,EPSONE,EPSCNEP, EPSCNPP , ERRORT , GAMBAR , CAMINF, 

E HCCMBTB ( 28 ) ,TTHC0M8( 7 ) , PHC OMB ( A) , NHCOMB , NPHCGMB ,HC T AB ( 2 8 I , 

G TTABHC ( 7 ) , PHC ( A ) ,NHC , NPHC , HET AB< 10 ) , TTABHE ( 10 ) , MFE , NHE , HWT AB ( 15 ) , 
I TTABHWl 15) , MHW,NHH, I ADJUST, I P LOT , L , MACHNO , M AX I TT , MDMA X , 

J MOOT 0( 20 ) , 

K MW02, MWSTR,NTP(7) ,PLTIME( 15) , PRAT (20) , PRFREQ, PSEXP , PS TAGTBI 1 0 ) , 

M TTPSTAG(IO) , MPSTAG.NPSTAG, PTMAX, PTMI N , GCT AB ( 1C I , TTABQC ( 10 ) , MOC , 

N NQC ,QRAT ( 20) , 

0 CRRAT (20 I »QRT AB (10) ,TTABQR( 10 ),MQR»NQR,R( 20) ,RIEXP,RNSI ,RO,RODP, 

0 R0P,RS(2C ) ,RSSM AX, S,STEBOL,T( 10,201 , T AUG , TBTAB ( 10 I ,TTABTB( 10 ) , 

5 MTR,NTB,TDPRIME,THETA(20> , TPR IME , XO, XORDE R ,ZS ( 2C ) , Z SMAX 
REAL MDCTC,MOOT,MCOOT, MSDOT, MV, STR ,MW02 , MACHNO, MOM AX 
DIMENSION AL(10) 

INTEGER S , SMI , SM2 

2H FORMAT (7E18.7) 

STATION ( 1,L ) X= L , ETA =0, 

W= 3.5* XC**2 * DELX I 
U= 8.3*H2(1,L)*H3<1,L) *X0 
XX = 3,0 * HI ( 1 , L) * H3 ( 1 , L ) * DELTA(L) 

SC= 3.0* CEL T A ( L ) **2 * DELETA 

G= -U*ROPCPP/W-XX*RGDPC/SC - H 1 ( 1 , L ) *H2 ( 1 , L I *H3 ( 1 , L ) *R 0*CP ( 1 , L ) 

1 /DELTAU 

PARTI = SC * OELETA 

E32L = (H1(2,L) + H1(1,L))*(H3(2,L)+H3(1,L) I * ( CKE T A ( 2 , L ) +C KETAC 1 , L ) ) / 
1(4,*(H2(2,L)+H2(1»L)H*9. 

E52L=(H1( 3,L )+Hl(2,L ) l*( H3( 3,L )fH3(2,L) ) * ( CKETA ( 3 , L ) +C KETAC 2 , L ) )/ 

1 (A.*(H2(3 ,L)+H2(2,L) ) ) 

D1LM32 = ( H2 ( 1 , L Ml ) +H2 ( 1 , LM2 ) ) *( H3 ( 1 ,LM1 ) +H3 ( 1.LM2 ) ) *( CKX I ( 1 ,LM1 )♦ 

1 CKX I (1.LM2 ) ) / (A.*(HI(I,LM1 )+H(l,LM2) ) ) 
D1LM12=(H2(1,LM1 )+h 2(1,LI)*(H3(L,LM1)+H3(1,L))*(CKXI(1,LM1 l+CKXI ( l 
l,L))/(A.*(Hl(l,LMl)+H1(l,L))) 

EPTA= ( — U*S I G P/ W -XX * SI GDP/ SC ) 

EPTB= EPT A * TB 
EPTA= EPTA* T ( 1 , L ) **3 
RSAVE = G 
C 

GO TO ( 1EG, 180, IROCOL 
15) CONTINUE 

BS1 ( 1 , L ) = -E32L/PART! 

C ( 1 , L ) = ( E52L + E32LI/PART1 
Z(L)= -E52L/PART1 
B ( 1 I = BS1 ( 1 , L ) + BSAVE + EPTA 
IF (IFIRST.EQ. 0 ) GO TO 180 


8A9COCCO 

85000000 

esioooco 

852CC0C0 
85300CC0 
B5AC0CCC 
85 5000C0 
85600000 
E570C0CO' 
85800001 
85900000 
86000000 
66100000 
86200000 
86300000 
86AOOOOO 
E650000G 
86600000 
6670C0C0 
86800000 
86900000 
87100001 
87100002 
87200000 
87300000 
67400000 
87500001 
87600000 
8770000G 
878C000C 
879C0QC0 
880C000C 
88100000 
882C00C0 
88300000 
88A000C0 
685COCOO 
886000C0 
887C00C0 
888000C0 
889300C0 
89000000 
89 ? 00000 
89200000 
893000C0 
89AC000C 
895000CC 
89600000 
89700000 
898C00C0 
89900000 
900C0CCO 
90 10C0C0 
902C0000 
90300000 
90 *-00000 
905C00CC 
90600000 
90700000 
90800000 
909C0C 00 
910C0C00 
91100000 
91 2C0CC0 
91 3CCQ0C 
91 ACCOOO 
9 15000 CO 
916C00C0 
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178 DC ( 1 ) = ( RSAVE - BS18(1,L))* T(l»L) -VB ( 1 I *T{). » L MZ ) - AB ( l , L I * 

3. T(liLMl) + EPTB 
GO TO 198 
180 CONTINUE 

VB ( 1 I =- D 1 LM32 / ( N*DcL X I ) 

AB(1,L)= (01LM32+ 9. G*DUM1 2 ) / { W*DELX I » 

3S1B(1,L)=-9.0*D1LM12/ <W*DELXI) 

IF ( IF IRST.EQ. 0 I GO TO 178 
19C B(L) = BSi B ( 1 , L ) + BSAVE +. EPT4 

DC ( L I = (BSAVE -BSIM.U >*T(1,L>- C < 1 , L I * T ( 2 , L ) -/(LI*T(3,L) 

1 +EPTB 
193 CONTINUE 

GC TO (202,30,: ) , IPOCOL 

STATION (M,L) X = L ETA GREATER THAN C, LESS THAN 1 

ENTR V COL XL M 
202 DO 210 M = 1 ,S 

2 L J ALIM) = F2(M,L)*H2(M,L)/H1.(M,L) 

W= 3.0 * XC * DELXI 
YY = DELTA(L) **2 * DELESQ 

DO 303 M= N1 » N2 
MM 1 = M-l 
MP1 = M+l 

XX = ETA(M» * ( A A ( L ) + A A( LM!) ) / ( 9-. 0* (DELTA { L I + DE LT A ( L Ml )) *DEL ET A ) 
XX1= ETAm * (AA(LM1)+AA(LM2I )/( 4. *( DE LTA < LM1 ) +DE LT A ( L 32 ) I * 
1DELETA) 

XY = 3.0* C(M,L) * HI ( M* L I / CKXHM.LI 
AN = £ TA ( M ) * A A ( L ) * CKXI ( M, L ) / DELT A ( L I 
AM = AN/ ( DELTA ( L I * DELESQ) 

G = H1(M, LI* H2 ( M , L I * H3(M,L) * RO * CP(M,L) 

AJ = AN / ( 4. 0* DELETA * XO * DELXI ) 

Ul= (H2(MP1,L)+H2(M*L) ) * ( H3 ( MP1 , L ) +H3 ( M , L ) > * ( ETA { MP1 ) +ETA( M ) ) 

1 /<4.3* (F1(MP1,L)+H1(M,L) I ) * A A( L ) 

U2= (H2(MM,L) + H 2 ( M , L ) ) * ( H3 ( MM1 , L ) +H3 ( M , L )) *( ETA( MM], ) +ETA ( M ) ) 

1/ (4.0* (HU MM 1 , L ) +H1(M,L) ) )*AA( L) 

DMLM32 = (H? (M.LM1 )+H2(M,LM2 > I *( H3 ( M ,L 31 ) +H3I M f LM2))*(CKXI(M,LMl) + 
1CKX I ( M , LM2 ) )•/ ( 4. * ( HI ( M,LM1 ) + H1 ( M , LM2) ) ) 

DMLM12=(H2(N» L.M1 )+H2(M»Lt)*(H3(M,LMi)+H3(M,LI)*(CKXI(M,LMl l+CKXI (M 
1»L) )/(4.*(Hl(M,LMl)+Hl(M,L) ) ) 

D3= -9.0*CMLM12* XX* ( T( MP 1 , L ) -T ( MM1 , L ) + 

1 T ( MP1 , LM3 ) - T(MM1,LM1) I 
D2 = 0MLM32 *(-XXl)* ( T( MP 1 , LM 1 )- T( MM1 , LM1 ) 

1 + T( MP1 » L M2 ) — T ( MM1 , LM2 ) ) 

DN =- (Di +D2I/W 

DN1 = AJ *( AL(MPl)* (3.C*T(MP1,L)-4.C*T(MP1,LM1)+T(MP1,LM2) ) 

1- ALIMM1 ) *( 3, 0*T ( MM 1, L ) - 4. 0*T( MM1 ,LM1 1 + T(MM1,LM2) ) ) 

BSAVE = -RCPCPP * XY/ W - G/OELTAU 
EPT4= -SIGP * XY /W 
EPTB= EPT4 * TB 
E P T4= EPT4 * T(M,L) **3 
C 

GO TO ( 240, 280 .IROOOL 
240 CONTINUE 

EMM12 = (FKM,L)*H1(MM1 ,L1) * ( H3 ( M , L ) +H3 < MM1 , L > ) * (CKETA ( M, L ) 

1 +CK.ET A( MM1 , L ) ) / (4.0* 1 H2 ( M , L ) +H2 ( MM 1 , L ))) 

E MP1 2= (Hi(M,L)+Hl(MPl,L) ) * ( H3 ( M , L ) +H3 ( MP 1 , L )) *( CKETA ( M,L ) 

1+ CKET A ( MP1 , L ) ) / (4.0* (H2(M,L) +H2 ( MP1 , L ) ) ) 

GH =G * ET A ( M ) * MOOT (L I / (CELT A { L ) *R0 *2.0* DELETA) 

A(M,L)= AM *U2 + EMM12/YY +GH 
C(M,L)= AM*U1 + EMP12/YY -GH 

BS1 (M,L)= AM* (-U1-U2 ) + (-EMM12 -EMP12J/YY 
B(M) = BSKM.L) + BSAVE + EPT4 
IF ( IF IRST.EQ. 0 ) GO TO 280 

278 DC ( M ) = DM + DN1 + BS AVE*T ( M , L ) - VB(M )*T(M,LM2) — AB ( M ,L ) * 

1 T(M,LM1»- BS1B(M,LI* T(M,LI + EPT8 
GO TO 300 


<91700000 
S18C0C3O 
91 9C30G0 
92000CCC 
92100C00 
92200000 
923000C0 
92400000 
92 SOOQCO 
92600000 
92700000 
928000C0 
9290000C 
93000000 
933 OOOCO 
9320OCCO 
93300 OCO 
9340C0CO 
93500CC0 
936C0C00 
937C0CC0 
93800CCC 
939C00C0 
940000CO 
941C00C0 
942COOCO 
94300000 
94400000 
945C0CC0 
946C00C0 
947C00C0 
94830000 
949C0CC0 
95000000 
95ICOOCO 
952C0000 
S53C0000 
S54C00C0 
955COOCO 
35600000 
95700000 
S58000CO 
959C0CC0 
96CC0C00 
96100000 
96200000 
96300000 
564CC000 
96500000 
96600000 
967C0000 
968000C0 
96900000 
97000CC0 
97100000 
97200000 
973C0000 
97400000 
975C0CC0 
97600000 
97700000 
57800CCC 
97900000 
98CC0CC0 
98100000 
982C00C0 
983C00CC 
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230 CONTINUE 

PART =W * XC * CELXI 
PART2=0MLM32/PART 
PART1= 9„C * 0MLM12/PART 
VB < M ) = - PART2 

AB(M,L > = PARTI + PART2 

BS1B(M,L) =- RART1 
IF (IFIRST.EQ.C ) GO TO 278 
290 B ( l I = 8S2B(M,L) + BSAVE + EPT4 

OC(L I = 0N + CN1f( BSAVE— BSi(M»L) ) * T(M,L)- A( M,L ) *T (MM1 ,L )- C(M,L) 

1 * T ( M PI , L ) + EPTB 
390 CONTINUE 

GO TO (3C1.80G I . IPOCCL 

STATION (S,L) XI =1, (X=L) , ETA=1, 

ENTRY C0LXL1 
331 CONTINUE 

W = 3. 3 * XCDXI 

WSQ = 3.0* XOXISQ 

DEDETA = CELT A ( L ) * OELETA 

TWDEL = 2.0* CELT A( L ) 

Ul = ( AACLM AAIL Ml ))/(2. *( DE LTA< L) + DELTA ( LMI ) ) ) 

U2=( AA(LM1 )+AA(LM2 I) / ( 2. * ( OE LT A( LMI I +DELT A < LM2 ) ) ) 

SP = ( HI ( S , L ) * XODXI+ 2. 0*TPRIME )/< HI ( S . L ) * XCOX I ) 

DHK = OELTA(L) * H2( S . L ) / CKETA(S,L) *SP 
DHK1 = DELTAILM! )* H2(S,LM1)/ CKETA(S.LMl) 

DHK2= DELTA (L M2 )* H2(S,LM2)/ CKETA(S,LM2) 

112 =3.3* CELE-TA * E(S,LI* DELTA I L ) * H2IS.LI* SP/CKET A( S.L I 
FF=1 • 3 / ( 3-. C*DEDETA**2) 

H = 8,0 * Hi ( S , L ) * D(S»L)/CKXI(S,LI 
PART = AAIL) /DEOETA 
AOD = PART/3.0 

A0D1 = (l.C + ETA(SMl)) * PART/2.0 

A0D2 = (HTA(SMl) + ETA ( SM2 ) ) *PART/2.C 

PART = 3.0* T(SM1,U-4.0*T(SM1,LM1)* TISM1.LM2) 

DSM32L = t H2 ( SM2 ,L)+H2(SMi,L I )*(H3< SM2 , L » +H3 I SMI ,L )) *( CK X I ( SM2 , L ) + 

1 CKX I ( SMI ,L) )/(4«*(Hl(SM2»L)+Hl(SMl,L) I ) 

PAR.T2=0SM3 2L*( 2. *T ( SM2 , L ) - 4. *T ( SM2 , L Ml ( + T ( SM2 , L M2 > + P AR T ) 

DSM12L= ( F2 ( SMI , L ) +H2 (S.L) ) * ( H3 ( SMI , L > +H3 ( S , L ) ) * (CKXKSMl.L) 

1 +CKXKS.L) )/(4.Q* (H1(SMI,L)+Hi( S.L) ) ) 

P A PT 1 = -9 i C*DSM12L*(3,9*T(S,L)-4. 0*T(S,LM1) + T(S,LM2)+PART) 

GSL = HUS. LI* H2 ( S ,L ) *H2 ( S , L ) * RO *CP(S,L) 

PARTW = -1.0/W + ACD r 

EPT4 = SI GP * H* P AR TW 

EPTB = EPT4 * TB 

EPT4 = EPT4 * T ( S t L ) **3 

DN = ACD * (PARTI + PART2)/(4.0* XODXI) + EPTB 
BSAVE=H*RC PC PP*( PARTW l-GSL/OELTAU 
GC TO ( 55'’, 650 ) , IPOCOL 
550 CONTINUE 

A J=GSL *MDCT ( L ) / ( R0*2 . 0*C EL TA ( L ) *QE LE T A ) 

DDSL= -FF*ZZ2 
QSAVE= DDSL* QS(L) 

ESM32L= (HI ( SM2 ,L I +H1 ( S Ml , L ) ) * (H3(SM2,L ) +H3 ( SMI , L ) ) * ( CKETA( SM2.L) 
1+ CKET A ( SMI «L ) ) / ( 4. C* ( H2( SM2 , L ) +H2 ( SMI , L > ) ) 

PARTE3=FF*ESM32L 
PARTD3 =. A0C*ADD2*DSM32L 
V ( L ) = -PARTC3- PARTE3- AJ 

ESM1 2L= (H1(SM1,L)+H1(S,L))*(H3(SMI,L)+H3(S,L) )*(CKETA(SM1 ,L) 

1 +CKETA(S,L))/(4.C* ( H2 ( SM 1 , L ) +H2 ( S , L ) ) ) 

PARTE1 = f F*9.C*ESM12L 
PARTD1 = ACC*ADC1*9.0*DSM12L 

A ( S , L ) = PARTDT + PARTD3 + PARTE3 + PARTE! + 4.0*AJ 
BS1 < S» L) = CPSL*SIG*T(S,L)**3 - PARTD1 - PARTE1 -3.C*AJ 
B(S)= BSKS.L) + BSAVE + EPT4 
IF (IFIRST.EQ.G ) GO TO 650 

648 DC ( S) = DN - VB(S ) *T(S,LM2) - AB < S , L ) *T ( S , LMI )- (BSIB(S.L) 

1 -BSAVE) * T(S,L)+ QSAVE + DDQSR 
GO TO 300 


9840000C 
9850000 0 
9 8 6 C 0 C 00 
98700000 
988C00CD 
589C0000 
99CC0CC0 
991 COC 00 
99200CC0 
993COOOO 
994COOOO 
99500000 
99600000 
997 C 0 C 00 
598000C0 
898CCCC0 
lOOCOCOOO 
1C01C0CC0 
1 002C 0000 
100300000 
100*00000 
100500000 
1006CO>000 
1 C07C000D 
1CC8C0C00 
1OO9COCO0 
lOlCOOCOO 
101100000 
1012C0GC0 
1D13C00C0 
101 400000 
IC15300C0 
i C160CC00 
1C17C0C00 
101 8C00C0 
101900000 
1 02000000 
1C21C0CG0 
102200000 
1 023 30 C 00 
102400C0C 
1C2500000 
102600000 
102700000 
1028COOCO 
102900000 
103000000 
103100000 
1032000C0 
1033C0000 
1034CCC0O 
1035000 00 
10 36 COC 00 
1 C 37C0C00 
l 03800000 
1039 300 CO 
1G4CCOOOC 
1041C0CCO 
104200000 
104300000 
1 C44Q0C00 
104500000 
1046C00CO 
104700000 
1048C OOOO 
1 049COOOO 
1050C0CC0 
105)00000 
1C5200000 
1C53C.0C00 
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650 CONTINUE 

UXODX I = W* X3CXI 

0SLM3=(H2(S t LMl M-H21 S,LM2) )*{H3(S«LMi) + H3(S»LM2 IS* (CKXI (St LMI ) 
HCKXIIS.L^I »/ <4„C* (Hl(S,LMl»+Hl<StLM2m 
DSLM1 = 9» 3* (H2 < S»L)+H2 <S,LM1 > ) *( H3 ( S » L ) +H3 I S , LMI ) ) *(CKXK$,LJ + 
1 CKXI (S.Lf'lI!/ <4,0*(Hl(S,LI+HMS,LMim 
QSLM1 = (-DSLM1 *U1 + DSLM3* U2» *DHK1/W 
QSLM2 * DSLM3*U2 *DHK2/W 
QSL=-U1* DFK* 0SLM1/W 

DDQSR= QSLM1* QSILMll +QSLM2 *QS(LM2) + QSL*QS(L) 
VB(S»=-OSLN3/WXeOXI+QSLM2*SIG*T(S,LM2>**3 
ABIS.L ) = {CSLM1+DSLM3)/WX0DXI +Q SLM1*SI G*T ( S , LMI ) **3 
BS1B(S,L)=-DSLN1/WX0DX I+QSL*SI G*T { S t L S **3 
IF (IFIRST.EQ.O ) GO TO 648 
69C B(L) = BSie(S,L> + BSAVE + EPT4 
DC(L l=ON+CSAVE * DOQSP 

1 - V( LI *T(SN2,L> - A(S,L) *TISMI,U - ( BS1 ( S , L ) -BSA VE (*T( S , L I 
800 RETURN 
END 


IC54C00C0 
1 C5500CC0 
1C560CC00 
1 C57000C0 

icsecooco 

105900000 
1 C6CC00C0 
1 06 1CC0C0 
1 C62C00C0 
106300000 
1C64C00C0 
1Q65CC000 
1 C66C0000 
1C6700000 
1C68QOOCO 
1C6S00P00 
IC70C00C0 
1C71C0C00 
1C72CC0C0 


Subroutine SQAERO.- Subroutine SQAERO computes convective and radiant heating 
rates and surface mass -loss rates and obtains variables which are functions of time, 
temperature, and pressure. The flow chart for subroutine SQAERO is as follows: 



39 


















noon on ooo 


The program listing for subroutine SQASRO is as follows 


SUBROUTINE SQAERO 

THIS ROUTINE COMPUTES THE HEATING RATES ANC THE MASS LOSS RATES 

COMMON /PICK/ A(1C,20),AA{20),A3{10,2C),ALPHA<20)»B(20{, 

2 BS1 1 1 0, 2C),BS1B( 10, 20, C 110,20) ,CB( 10,20) ,CK( 10) ,CKET At 1C, 20 » , 

A CKXI (10,201 , COST (20 1 ,CP ( 10, 20 5 , D ( 10 , 20 1 , OC ( 20 ) , 

6 OELESQ,OELETA,DELTA<20) ,DELX I , DEL XI SG , E (1 0 ,20 5 , E IGHT3 » 

8 ELAM( 20) , ETA( 10 1 ,EXPG,F (10,20 I , GG, G IMACH , HI ( 10 ,20 ,F2(1C,20) , 

A H3 ( 10 , 20 ) » HC( 201 ,HCOMB< 2C 1 , HE , HW I 20) , I F IRS T, I POCCL , ITC , ITR , ITT, 

C ITTO,LMl,LM2,MCDOT(2C),MDOT(221,MSOOT(2C! , PI 02, PR ELCC ( 20) ,QC(20) , 
E 0C1, QC.NET (20) ,QCCM8(20) ,QR(20 ) ,QR1 ,QS(20) , FNS , RCDPC ,ROPCPP, 

G RSS ( 22 I ,RST02 ,SIG,SIGOP, SIGMA, SI GP, SINT (20) , SM1,SM2 ,TAU,TB* 

1 TT(10,20) ,TTF (10,20) ,TWDELXI,TWQGI»V(20) *V8(10) ,X(22! ,XOXISQ, 

K XCDXI,Y(iO,20),Z(20) ,ZB(10) 

COMMON /INPUTS/ DUMMY (I) , AEXP, ALCTABI 10 ) , TT ALC < 10 ) , MAC PHC,NALPHC , 

2 ALPHA T( 10) ,TALPHA(1C) ,MAL PH A , NAL PHA , AL ST AB ( 10 ) , TT AL S ( 10 ) , MALPHS , 

4 NAL PH S, AS EXP, BETA, BEXP, BSEXP , CE »CKET ATB ( 50 ) , T 7CKETA ( 10) , 

6 ETATA8I5) , NCKETA , NETA ,CKX I TAB ( 50 ) , TTCKX I (10 ) , X I TAB ( 5 ) ,NCKXI , 

8 NX I , CORDS Y ,CPCP,CPP,CPTAB(10) ,TTABCP( 10) , MCP ,NCP ,DELTA0(20) , 

A DELTAU,DELTMIN,0TMAX,ELAMTB(28) , TTELAM( 7 ) , PEL AM( 4 ) , NELAM, 

C NPELAM,ENDTAU,EPS0NE,EPS0NEP,EPS0NPP,ERR0RT,GAM8AR,GAMINF, 

E HC0M3T8I28) , TTHCCMB ( 7 ) , PHCCMB ( 4 ) , NHCCMB , NPHCCMB ,HCT AB (28 I , 

G TTABHC(7),PHC(4) ,NHC ,NPHC,HETA3( 10) , TTABHE ( 10 ) , MHE , NHE »HWT AB ( 15 ) , 
I TTABHWU5) , MHW,NHW, I ADJUST, I PLOT , L , M ACHNO , MAX ITT , MDMAX , 

J MDOTO ( 20 ) , 

K MWC2, MWSTR,NTP(7) ,PLT IME( 15 ) , PR AT ( 20 ) , PRFRE Q , P SE XP , PS TAGTB( 1 0 ) , 

M TTPST AG( 10) ,MPSTAG,NPSTAG,PTMAX,PTMIN,QCTAB(10) , T TABQC ( 10 ) , MQC, 

N NGC , QRAT (20 ) , 

0 QPRAT ( 20 ) ,QRT AB (10) , TTA60R (101 , MQR , NCR , R ( 20 ) , R I EXP , RNS I ,R0, ROCP , 

Q RCP.RS12C ),RSSMAX,S,STEBCL,T( 10,20) ,TAUO*TBTAB( 10 ) , TTABTB( 10 ) , 

5 MTB,NTB,TDPRIME,THETA(20) , TPRI ME , XO, XCRDER , ZS ( 2C ) , ZSMAX 
REAL MDOTC *MDOT , MCDOT , MSDOT , MW STR , MW02 , MACHNO , MDMAX 
INTEGER S,$M1,SM2 

LOOK UP CP, CPBAR, CKN ,ETC. AS FUNCTIONS OF TEMPERATURE 
DO 11 N=1,L 

CALL FTLUP (TT (S,N), ALPHA(N) ,M ALPHA, NAL PHA , TALPHA , ALPH AT ) 

11 CALL FTLUP ( TT ( S ,N ) , HW ( N ) , MHW , NHW , TT ABHW , HWT AB ) 

IF (ITT.NE.l) GO TO ICO 
LOOK UP FUNC TICKS OF TIME 

CALL FTLUP (TAU,ALPHAC,MALPHC,NALPHC,TTALC,ALCTAB) 

CALL FTLUP ( T AU , ALPH AS , MAL PHS, NAL PHS ,TT AL S , AL STAB I 
CALL FTLUP (TAU, HE, MFE, NHE, TTABHE, PETAB) 

CALL FTLUP (TAU,PSTAG,MPSTAG,NPSTAG,TTPSTAG,PSTAGTB) 

CALL FTLUP (TAU.QCl , MQC, NQC, TT ABQC , QCT AB ) 

CALL FTLUP ( TAU , QR 1 , MCR , NOR , TTA8QR ,QRTAB> 

CALL FTLUP { T AU , TB ,MTB ,NTB , T T ABTB , TBT AB I 
TB =T3**4 

ADJUST CONVECTIVE AND RADIANT HEATING RATES AND THE PRESSURE AND 
HEATINGOISTR I BUT I CN TO SHAPE CHANGE (ADJUST QC1 ,QR1 , PRAT, QRAT ) 

IF (CORDSY.NE.O ) GO TO 20 
CALL ADJUST 
20 00 30 N=1,L 

DELTAO(N) =DELTA(N! 

QR ( N ) = QR 1 * QRRAT(N) 

QC ( M = QC1 *QR AT (N I 
PRELOC(N) = PSTAG 4 PRAT(N) 

CALL D ISCCT ( TT ( S ,N ) , PRELOCI N ) .TTABHC , HCTAB , PHC , 11 , 28 , 4,HC ( N) ) 

CALL DISCCT(TT(S ,N) , PRELOCI N) ,TTELAM , EL AMTB , PEL AM, 1 1 ,28, 4, EL AM (N ) 

1 ) 

31 CALL DISCCT ( TT ( S ,M .PRELOC ( M , TTHCCMB, HCOMBTB , PHCCMB, 1 1 , 28, 4, 

1 HC0M8 (N) > 


1C 7300000 
107400000 
1C7500000 
107600000 
1 C77C0000 
107800000 
1C79C0000 
1C8C00001 
IC8100000 
1082C0000 
108300000 
1C8400000 
1C8500000 
108600000 
1087C0000 
108800000 
1C89C0000 
109000000 
1C9100000 
109300001 
1C9300002 
109400000 
109500000 
109600000 
1C97C0001 
109800000 
1C9900000 
110000000 
110100000 
1 10200000 
110300000 
110400000 
11C500000 
1 106C0000 
1 107000C0 
1108C00C0 
110900000 
1110C0000 
111100000 
111200000 
111300000 
111400000 
1115C0C00 
111600000 
1 11700000 
111800000 
111900000 
112000000 
112100000 
U22C000C 
1 12300000 
l 124C0000 
1125COOOO 
112600000 
112700000 
1 128COOCO 
112900000 
113000000 
1131COOOO 
113200000 
1 133CCOOO 
1 134000C0 
1 135000C0 
113600000 
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C COMPUTE QS ACROSS FRCNT SURFACE 
BAT = L,0 - BETA 
133 DO 200 N=T ,L 

CELL = HE /GC ( N ) 

CAT = QC(N) * (l.C - HW(N)/HE) 

BL0CK= ( ALPHAC *MCDOT(NI + ALPHAS *MSDGT ( N ) ) * CELL 

QCNET(N) '= CAT *(1.0 - BAT *(0.6* BLOCK - C.084 * BL0CK**2) 

1- BETA * BLOCK I 
QCCMB(N)= MCDOTt N 1 * HCOMB(N) 

QS ( N ) = QCNET(N) + ALPHA* OP(N»- MSDOT (N i *HC ( N ) + CCOMB(N) 

2)0 CONTINUE 
RETURN 
C 

C THIS PART OF ROUTINE CCMFUTES ROOTS 
C 

ENTRY SQAEPOM 
DO 1030 N=1,L 
C 

C COMPUTE MSDOT MASS LOSS RATE DUE TC SUBLIMATION 

C 

IF (ASEXP ) 310,3X5,313 
335 MSDOT ( N I =C .0 
GO TO 330 

310 BLOCK =-BSEXP/TTF(S,Nl 

MSDOT ( N) = ASEXP * PREL0C1NI **PSEXP * EXP ( BLOCK )*R ( 1 ) **RIE XP 
333 COLL = (HE-HW(N))/(QCNET(N)*ELAM(N)) 

C 

C COMPUTE MCDOT MASS LOSS RATE DUE TO OXIDATION 

C 

C HALF ORDER OXIDATICN 
C 

333 IF (AEXP) 390,385,390 
385 MCDOTIN) =0.0 
GO TO 330 

333 MCDOTIN) = AEXP * EXP ( -BEX P/TTF( S ,N 1 ) 

IF (X0RDER-C.5)' 9CC ,403,600 
433 ABC = 4.3* MCDCT ( N I **2 * PRELOC(N) * CE * RST02 ' 

PART = COLL * MCDOT ( N ) **2 * PRELOC(N) * RSTC2 

TEST = ABC/ PART **2 

IF (TEST.LT.7.E-12IGO TO 420 

MCDOTIN) =• 5*( (-PART ) + SQRT (PART**2 + ABC)) 

GO TO 900 

420 MCDOT(N) = CE /COLL 
GO TO 930 
C 

C FIRST CRDER OXIDATION 
C 

630 MCDOTIN) = MCDOT ( N ) * PRELOC(N)* RST02 * CE/I1.0 + MCDOTJ N)*PRELOC 
1 (N)* C0LL*RST02 ) 

C 

C MOOT IS EQUAL TC TFE LARGER OF MSOCT AND MCDOT 
C 

930 IF (MCDOTIN). LT. MSDOT(N) ) GO TO 950 
MOOT ( N I = MCDOTIN) 

MSDCT ( N ) = 0.0 
GO TO 100C 

950 MOOT ( N )= MSDOT (N I 
MCDOTt N) = C.C 
10)0 CONTINUE 
RETURN 
END 


1 137C0000 
1 138000C0 
113900000 
114000000 
114100000 
114200000 
114330000 
1144X0000 
1 145C0000 
1 146C0000 
1 14700000 
1 14800000 
114900000 
1150C00C0 
115100000 
115200000 
1 153C00C0 
115400000 
1 15500000 
115600000 
115700000 
1158C0000 
115900000 
116000000 
1161C00Q0 
1 16230000 
116300000 
116400000 
116500000 
116630000 
1 167C00C0 
116800000 
116900000 
1170C0C00 
117100000 
117200000 
117300000 
117400000 
117500000 
1 1 7600000 
117700000 
1 17800000 
117900000 
118000000 
118100000 
118200000 
118300000 
118400000 
118500000 
118600000 
118700000 
118800000 
118900000 
119000000 
119100000 
1 19200000 
119300000 
119400000 
1 195C0000 
1 196C0000 
119730000 
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he 

The flow chart for subroutine ADJUST is as follows: 












oooo O uou 


The program listing for subroutine ADJUST is as follows: 


SUBROUTINE ADJUST 

THIS ROUTINE ADJUSTS THE CONVECTIVE AND RADIANT HEATING RATES, THE PRESSURE 
AND HEATING DISTRIBUTION TO SHAPE CHANGE (ADJUST «C1 ,QRl ,PRAT,QR AT ) 

COMMON / PICK/ A(10,20) , A A ( 20 > , A8 ( 10 ,20 ) , ALPH A( 20 ) , B I 20 > , 

2 BS1(10,2C),BS1B(1G,20»,C(10,20).CB(10,20),CM10>,CKETA( 10,20) , 

4 CKXIJ 10,20) ,COST( 20! ,CP( 10, 20) ,0(10,20) ,DC( 20) , 

6 DELESQ,DELETA,DELTA(20),DELXI,DELXISQ,E(10,20) ,EIGHT3 , 

8 ELAM( 20) , ETA ( 10) ,EXPG,F( 10,20 ) ,GG,GIMACH,H1 ( 10,20 J,F2 (10,20) , 

A H3(10,20)»HC(20)»HCOMB(2C),HE»HW(2C),IFIRST,l ROCOL ,ITC,ITR,ITT, 

C ITTG,LMl f LM2,MCD0T(20),MD0T(22l, MSDOT (20),PID2,PREL0C(20),QC(20), 

E QC1 ,QCNET ( 20) ,QCOMB(20) ,QR ( 20 ) ,QR1 ,QS( 20) , RNS ,RCDPC ,R0PCPP » 

G RSS122) ,RST02,SIG,SIGDP,SIGMA,SIGP,SINT{20),SM1,SM2,TAU,TB, 

1 TT(10,20),TTF (10,20) , TWDELX I , TWOGI , V( 20 ) , VB( 10 ) ,X (2 2 ) ,XDX ISQ, 

K XCDXI,Y(10,20),Z(20),ZB(10) 

COMMON /INPUTS/ DUMMY ( 1 ) , AEXP, ALCT AB ( 10 ) , TTALC ( 10 ) ,MAL PHC, NAL PHC, 

2 ALPHAT(IC) ,T ALPHA (10) ,MAL PHA, NALPHA, AL STAB ( 10 ) ,TTALS( 10) .MALPHS, 

4 NALPHS,ASEXP,BETA,BEXP,BSEXP,CE,CKETATB(50),TTCKETA(10) , 

6 ETATA8(5 ) ,NCKETA, NETA ,CKX I TAB ( 50 ) ,TTCKXI( 10 ) ,XITAB( 5) ,NCKXI , 

8 NXI ,C0RDSY,CPDP,CPP,CPTAB(10) , TT ABCP ( 10 ) , MC P , NCP, CELT AO ( 20 ) , 

A DELTAU,DELTMIN»DTMAX»ELAMTB(28I , TTELAM ( 7 1 , PEL AM (4 ) , NE LAM, 

C NPELAM,ENDTAU,EPSCNE, EPSONEP , EPSONPP,E RRORT , GAMBAR , GAM INF, 

E HCOMBTB (28) , TTHCCMB ( 7 ) , PHCOMB ( 4) , NHCOMB , NPHCOMB ,HCT AB ( 2 8) , 

G TTABHC(7) ,PHC(4) ,NHC, NPHC , HET AB( 10 I ,TT ABHE ( 10 > , MFE ,NHE, HWTAB ( 15), 

I TTABHW( 15 ) ,MHW,NHW, I AD JUST , I PLOT ,L , MACHNO , MAX ITT , MDM AX , 

J M00T0I20 ) , 

K MW02,MWSTR,NTP(71 , PLTIME (15 ) , PRAT (20) , PRFREQ , PSEXP , PST AGTB( 10 I , 

M TTPSTAG(10),MPSTAG,NPSTAG,PTMAX,PTMIN,QCTAB( 10),TTAB0C( 10) ,MQC, 

N NQC ,QRAT ( 20.) , 

0 QRR AT ( 20 ) ,QRT AB ( 10 1 , TTABQR( 10 ) , MQR ,NQR »R( 20 I ,R I EXP ,RNSI ,RO,RCDP, 

Q ROP ,RS( 2C ) ,RSSMAX,S,STEBOL,T( 10,20) ,TAUC ,TBTAB( 10 ) , TT ABTB( 10 ) , 

5 MTB,NTB,TDPRIME»THETA(20)»TPRIME»X0,XQRDER,ZS(20),ZSMAX 
REAL MD0TC,MD0T,MCD0T,MSD0T,MWSTR,MHQ2 ,MACHNO,MDMAX 
INTEGER S,SM1,SM2 

DIMENSION PS I ( 20 ) 

DIMENSION UEUK20), AL( 20 ) , A I NT ( 20 ) , YY ( 3 ) 

NS PI = NSTEP + 1 
’DO 50 N=1,L 

RSS(N) = RS ( N ) + DELTA(N) *CCST(NI 
50 ZS ( N ) = ZS ( N ) + (DELTAO(N) - DELT AIN))* SINT(N) 

IF ( IADJUST.EQ.O ) RETURN 

RNS= ( ZS ( 2 ) **2 * R S S ( 2 ) ** 2 -2. 0*ZS ( 2 ) *ZS ( 1 ) + ZS(1)**2)/ 
1(2.0*(ZS(2)-ZS(1))) 

SQRNS = SORT ( RNS ) 

ADJUST RATE TO SHAPE CHANGE 

QC1 = QC1 * SORT ( RNSI/RNS ) 

QR1 = QR1 * RNS/ RNSI 
PSI (11=0. 

M=1 

100 DO 200 N=2,L 
NP1 = N+l 
NM1 = N-l 

IF (N.EQ. L) GO TO 130 

TANPHI = (RSS(NPl) -R SSI NM1 ) ) / ( ZS ( NP 1 )- ZSINM1H 
GO TO 150 

130 T ANPHI = (RSS(L l-RSSUMl ) ) /( ZS(L )-ZS(LMl ! ) 

150 PHI =• ATAN (TANPHI ) 

PSI ( N) =P I D2-PH I 
200 CONTINUE 

NEW PRESSURE DISTRIBUTION 
DO 250 N=1,L 

PRAT(N) =(1.0 - GTMACH) *COS( PSI ( N! 1**2 + GIMACH 

UEUI(N) = SORT ((1.0+ TWOGI) *< 1. O-PRAT ( N ) **EXPG) ) 

250 CONTINUE 

OBTAIN NEW FEAT DISTRIBUTION 


119800000 
119900000 
120000000 
120100000 
120200000 
120300000 
120400000 
120500000 
120600001 
120700000 
120800000 
120900000 
121000000 
121100000 
121200000 
121300000 
121400000 
1 21 5C0000 
121600000 
121700000 
121 9000C1 
121900002 
122000000 
122100000 
1222C0000 
122300001 
122400000 
122500000 
122600000 
122700000 
122800000 
122900000 
123000000 
123100C00 
1 23200000 
123300000 
1234C0000 
123500000 
123600000 
123700000 
123800000 
1 239C0000 
124000000 
1 241 COOOO 
1242C0C0C 
124300000 
124400000 
124500000 
1 246GC000 
1 247COOOO 
124800000 
1249C0000 
1 25C00000 
125100000 
125200000 
1253COOCO 
125400000 
1 255C0000 
125600000 
125700000 
125800000 
125900000 
126000000 
126100000 
126200000 
126300000 
1264COOOO 
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C EVALUATE INTEGRAL AT L =0 
AL ( 1 ) =0* 0 

AINT0=PRAT( 1 )*UEUI( 1 )* RSS(1)**2 
Z 7 0 CGNTINUE 

OR AT ( 1 1 = 1® 0 
DC 600 N = 2 , L 
NM 1=N- 1 
NP2 =N-2 
A I NT = A I NTC 
SUMH1=0» 

IF ( N, EQ« Z ) GO TO 310 
DO 300 1=2, NM1 
300 SUMHI=SUMH1+H1 (S,I ) 

310 AL ( N )= X t 2 > *1 SUMH1 + »H1IS,1)+ H 1 ( S, N) ) /2. 0) 

EVALUATE INTEGRAL 

IF CN.EQ. 2) GO TO 500 
C EVALUATE Y ( 0 ) , Y (1 ) , Y ( 3 I 
DO 400 K= l.,3 
NMK = N- ( 3— K » 

400 YY(K)= PRAT(NMK)*UEUI ( NMK )*( RSS(NMK)*#21 
C0EF2= ALINM2)- ALIN) 

POXO= I AL ( NM2 ) - AL(NMl))* CCEF2 

PI Xl= I ALI NM1 ) - ALINM2I)* (AL(NMl)- ALIN)) 

P2X2= ( AL ( N ) - ALINM2)) * ( AL ( N)-AL ( NM 1 ) ) 

COEF 1= <3.0* ALINMD-2.0* ALINM2) - AHN)»/POXC 
CCEF3=(2.C*AL(M + ALINM21- 3,0* ALINMl))/ P2X2 

AINTIN) = ( I AL ( N ) - AL ( NM2 ) ) **2 /6, 0 ) * ( YY(1)*C0EF1 + YY(2)*C0EF 2/ 

1 P1X1 + YY 1 3 ) * C0EF3 ) 

IF (N.GT.2) AINT (N) = AINT I NM2 ) + AINTIN) 

GO TO 590 

C N= 2 

530 YY 1 2 )= I PRAT (1 ) + PRATI 2) >* I UEU 1 1 1 ) + UEUII2) ) * f (RSSI1I+ RSSl2))/2.0 
1 ) **2 /4,G 

YY I 3 )= PR ATI 2 ) * UEUH2) *(RSSI2)**2) 

AINT(N)=AL(2)*(4.C* YYI2) + YY(3))/6,C 
590 ANUM=PRAT<N)*UEUKN)*RSS<N) *SQRNS 
QRATIN) = ANUM / I SORT I A INT I N ) ) * GG ) 

600 CONTINUE 
RETURN 
END 


126500000 
1266000C0 
126700000 
126800C00 
126900000 
127000000 
127100000 
127200000 
127300000 
127400000 
I 275C0000 
127600000 
1277C0000 
1 278000C0 
127900000 
128000000 
128100000 
1282 00000 
1283C0000 
128400000 
128500000 
128600000 
128700000 
128800000 
128900000 
129000000 
129100000 
129200000 
129300000 
129400000 
1295C00C0 
129600000 
1297C0Q00 
129 800000 
1295C00C0 
130000000 
1 30100000 
130200000 
130300000 
130*00000 
130500000 
130600000 


Subroutine Z PRINT . - Subroutine Z PRINT writes the output data. The flow chart for 
subroutine Z PRINT is as follows: 
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The program listing for subroutine Z PRINT is as follows 


SUBROUTINE ZPRINT 

COMMON /PICK/ 6(10,20) ,AA( 20 ) , AB ( 10 , 201 , ALPH A! 20 ) ,B( 20 ) , 

2 BS1 (10,20),BS1B(10,20)»C(10,2Q),CB(10,20)»CK(10)*CKETA(10»20I , 

A CKXI( 10,20) ,COST( 20) , CP (10,20 ),D( 10,20 ) ,OC (20) , 

6 DELES Q,DELETA, DELTA (20) , OELX I , DELX ISQ , E ( 10,20) , E IGHT3 , 

8 EL AM{ 2 J ) , ETA ( 10 ) , EXPG ,F ( 1 0 , 20) , GG.GIMACH ,H1 ( 10 ,2C ) , H2 ( 10, 20 ) , 

A H3(10,20),HC(20), HCQMB( 2C I , HE , HW ( 20 ! , IF IRST , IROCOL , I TC , I TR , I TT, 

C ITT0,LM1,LM2,MCD0T(20)»MC0T(22)» MS DOT ( 20 > ,P ID2 , PPFLGC I 20 ) ,0C ( 20 ) , 
E QC1 »QCM£T (20 ) ,QC0MB(20S ,QR ( 20 ) , QR1 ,QS ( 20 ) ,RNS,RODPC ,ROPCPP, 

G RSS ( 22 ) ,RSTO2,SIG,SIGDP,SIGMA,SIGP,SINT(20),SMl,SM2,TAU,TB, 

1 TT(10,20)»TTF(10,20)»TWDELXI*TWOGI*V(20) , VB( 10) , X( 22 ) ,XDXI SQ , 

K XODXI ,Y(10,20) ,Z(20J ,ZB(10) 

COMMON /INPUTS/ DUMMY! 1) , AEXP , ALC TAB ( 10 ) ,TTALC(10) , MAL PHC, NAL PHC , 

2 ALPHAT ( 10 ) ,T ALPHA (10) ,MAL PH A, NALPHA, ALST AB ( 10 ) , TTALSC1C ) , MALPHS , 

4 NALPHS, ASEXP, BETA ,BEXP, BSE XP ,CE,CKETATB (50 ! , T TCKET A ( 10 ) , 

6 ETATABO) ,NCKETA,NETA,CKXITAB(50 ) ,TTCKXI( 10 ) ,XITAB( 51 ,NCKXI , 

8 NXI ,C0RDSY,CPDP,CPP»CPTAB(10) , TT ABCPUO ) , MCP , NCP .DELTAC ( 20 ) , 

A CELTAU»DELTMIN»DTMAX,ELAMTB(281» TTEL AM ( 7 ) , PELAM( 4 ) , NELAM, 

C NP£LAM,ENDTAU,EPSCNE,EPSONEP, EPSONPP, ERRORT, GAMBAR , 6AMINF, 

E HC0MBTB(28) ,TTHCCMB(7 ),PHCGMB( 4) .NHCOMB , NPHCOMB ,HCT AB ( 2 8 ) , 

G TTABHC(7 I ,PHC(4) ,NHC , NPHC , HET AB( 13 I , TT ABHE ( 10 ) , MFE, NHE, HWTAB ( 15), 
I TTABHH(15) ,MHW,NHW, I ADJUST, I PLOT , L , MACHNO , M AX ITT , MDM AX, 

J MDOTQ ( 20 ) , 

K MW02,MWSTR,NTP(7) , PLT IMF < 15 ) , PR A T ( 20 ) , PRFRE Q , PSEXP , PSTAGTBI 1 0 ) , 

M TTPSTAG(10),MPSTAG,NPSTAG,PTMAX,PTMIN,QCTAB(1C) »TTARQC( 10) ,MQC, 

N NOC ,QRAT ( 20 ) , 

0 CRRAT (20 I ,QRT AB(10),TTAB0R( 10 ),MQR,NQR,R( 20) ,PIEXP,RNSI ,RO,ROOP, 

Q ROP ,R S (20 »R3SMAX,S,STEBCL,T(10,2Q)»TAUG,TBTAB(1C ) , TT ABTB( 10 ) , 

5 MTB,NT3, TOP RIME, THETAI20) ,TPRIME ,X0,X0RDER,ZS(2C) ,ZSMAX 
REAL MOCTC,MDOT,MCDOT,MSDOT, MW STR , MW02 , MACHNO, MDMAX 
INTEGER S , SMI , SM2 

DIMENSION QRRI20) 

EQUIVALENCE (QRR( 1) ,HH 1, 1 ) ) 

‘DO 10 N=1 , L 

10 QRR ( N ) = SIG * TTF(S,N)**4 
WRITE (5, 98) 

98 FORMAT ( *0*) 

WRITE (5,100) TAU.DELTAU 
110 FORMAT <*CTAU=*F1C.4,14X*DELTAU=*F9.6) 

WRITE ( 5 , 1C1 ) QC1, QR 1 , HE 

101 FORMAT (*C*14X*QC=*E11.4,5X,*QR=*E11.4,5X,*HE=*E11.4 ) 

WRITE (6,1021 T ( S , 1 ) 

102 FORMAT (15X*T( S,1)=*E11.4) 

WRITE (6,105) 

105 FORMAT ( * C*14X* TEMPER A TURE ( M ,N) * ) 

WRITE (6,nC)(X ( N ) , N= 1 , L I 
110 FORMAT (* ETA*6X*X=*15F8.5/(12X,15F8.5( ) 

DO 115 M=1,S 
MM= S- (M-l) 

115 WRITE (6,120) E TA( MM ) , ( TTF ( MM ,N ) ,N=1 , L ) 

120 FORMAT ( F6. 3 ,6 XI 5F8, 1/ ( 1 2X , 1 5F8. 1 ) 1 

140 FORMAT (* ETA*6X*X=*1C(F9. 5,3X»/( 12X,10(F9. 5,3X) ) ) 

150 FORMAT ( Ft, 3 , 6 X10E12 . 4 / ( 12 X , 10 El 2 , 4 I I 

WRITE (5,155) 

155 FORMAT ( *C*l4X*MD0T( Nl — SURFACE MASS LCSS RATE*) 

WRITE (6,140) (X (Nl,N=l,L> 

WRITE (6,150) ETA(S) ,(MDCT(NI ,N=1 ,L) 

WRITE (6,165) 

165 FORMAT ( *C*14X*MCD0T( N ) — SURFACE MASS LOSS RATE DUE TO OXIDATION*) 
WRITE (6,150) ETA( S » , (MCDOT (N ) ,N = 1 ,L ) 

WRITE (6,170) 

1 7 0 FORMAT (*C*14X*DELTA(NI — MATERIAL THICKNESS*) 

WRITE (6,150) ETA(S), ( DELTA! N ), N=1 ,L ) 


130700000 
1 308C00C0 
130900000 
131000000 
131100000 
131200001 
1313C0000 
131400000 
131500000 
131600000 
131700000 
131800000 
131900000 
132000000 
132100000 
132200000 
132300000 
132500001 
132500002 
132600000 
132700000 
1328C0000 
132900001 
133000000 
133100000 
133200000 
133300000 
1334C0000 
133500000 
133600000 
I 33700CC0 
1338300C0 
133900000 
134000000 
134I000C0 
1342C0000 
134300000 
1 344000C0 
1 34500000 
134600000 
134700000 
1 34800CC0 
134900000 
1 35000000 
135100000 
1352C0000 
1353C0000 
135400000 
135500000 
1 35600C00 
1 35 7COOCO 
135800000 
1 35900CC0 
1 36000000 
136100000 
1 362000G0 
1 363C0000 
136400000 
136500000 
1 366C0000 
1367C0000 
1368CC0C0 
136900000 
1 37COCOOO 
1 37100000 
137200000 
137300000 
1 374C0000 
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WRITE (6,175) 

175 FORMAT (*C*14X*QRAT(N) — -RATIO OF LOCAL HEATING TO STAGNATION HEATI 
1NG* ) 

WRITE ( 6 50 ) ETA(S),( QP AT ( N) »N=1 , L ) 

WRITE! 6,176) 

176 FORMAT(*0*14X*PP ATINl — RATIO OF LOCAL PRESS TO STAG PRESS* I 
WRITE (6,150 ETA(S) , (PRAT(N) ,N=1,L) 

WRITE (6,180) 

130 FORMAT (*C*14X*GS(N) — NET HEAT INPUT*) 

WRITE (6,150) ETA(S), (QS (N|,N=1,L) 

WRITE (6,190) 

190 FORMAT ( *C*14X*QRR (N ) — RER AO I ATION* ) 

WRITE (6,150) ETA(S), (QRR(N) ,N=1,L) 

WRITE (6,200) 

200 FORMAT (*C*14X*QCCM8 ( N ) — HEAT DUE TO COMBUSTION FOR OXIDATION*) 
WRITE (6,150) ETA(S), (QCOMB ( N ) ,N=1 , L ) 

WRITE (6,400) I TC , ITR , ITTO .IROCOL 

400 FORMAT ( *C NO. ITER. C0L. = *I4 , 5X, *N0. ITER. ROW = * I 4, 5X»*T0TAL 

1 NO. ITER. = *I8, 5X*IR0CCL=*I 3) 

RETURN 

END 


1 3 75000C0 
137600000 
1 377C000C 
13780000C 
1 37500000 
1 380QGOCO 
138100000 
1 3823000C 
138300000 
138400000 
138500000 
1 38600CCC 
138700000 
138800000 
138900000 
139000000 
139130000 
1 39200000 
139300000 
1 394000C0 
139500000 
1 39600000 
1 397C0000 
139800000 
1 399000C0 
140000000 
140100000 


Subroutine SOLMAT .- Subroutine SOLMAT solves a system of linear equations in 
which the matrix of coefficients is a tridiagonal matrix. The method of solution is equiv- 
alent to Gaussian elimination. The flow chart for subroutine SOLMAT is as follows: 
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The program listing for subroutine SOLMAT is as follows: 


SUBROUTINE SOLKAT< A,B,C» Z, V, D,T »N! 

01 HENS ION Ml 20! ,SVJ20S S GC20 8 9 TI20 8 * A C20 I S B 120 ! , Cl 20 8 f D { 20 5 

COMMON /HCLO/ THIN 

C 

C THIS ROUTINE SOLVES THE TR IDI AGONAL {EXCEPT TMO ELEMENTS 5 MATRIX 

C 

M { 1 !=B { 1 1 

svm= cm i bi i ! 

X= Z/BJll 
GUS® DUI/Wll ! 

NMi=N“l 
NM2 = N— 2 

00 100 K=2 S N 
KM1 = K-l 

IF ( K* EQ® M GO TO 20 
H«K8 = B«tO - A{ K 5 *SV { KM1 } 

IF IK«EQ.28 GO TO 10 

4 SV (K 8= CJK8/WCKS 

5 G(K) = SD5KI- A<KJ*G(KM188/WIK8 
GO TO 100 

10 SV(2 8 =(C(2)-X*A{2)8/W{28 
GO .TO 5 

20 M(N)» B( N 1- < A { N 8- V*SV(NM28 >*SV< NM18 

30 G(N8=<D(N)-A<N)*GIKM18-V*G<NM28*V*$VINM28*G(KM18 J/W{ N8 
100 CONTINUE 
TCN)*GCN8 
DO 200 K=1,NM2 
KK= N-K 

T(KK»= G{ KK 8— SV { KKI *T (KK+ 1 8 
200 CONTINUE 

TC 1 8 = GUI- SV(lt*T(28- X*TI38 
IF (TMIN® EC®0« 1 RETURN 
DO 300 1 = 1 S N 

1 FIT ( I 8 .LT.TMIN8 T(I8=TMIN 
300 CONTINUE 

RETURN 

END 


PROGRAM INPUT, OUTPUT, AND DIAGNOSTICS 


Input 

Examples of input data are given in appendix B. The first card of the input is iden- 
tification for the job. Any identification may be written in column 1 to and including 
column 80. 

FORTRAN IV NAMELIST -with the name D2430 is used to load the input data. Each 
input variable is initially set equal to zero by the program unless otherwise stated. 

At least four inputs are associated with each table input: the dependent -table values, 
the independent -table values, the number of entries in the table, and the order of interpo- 
lation. The number of entries in the dependent and independent table must be the same. 
This is specified by a FORTRAN variable beginning with the letter N. The order of inter- 
polation is a FORTRAN variable beginning with the letter M and may be 0, 1, or 2. For 
example, for first-order interpolation of the specific -heat array, set MCP=1; for second - 
order interpolation, set MCP=2. If the specific heat is a constant, set MCP=0. 
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1*80200000 

140300000 

140400000 

140500000 

140600000 

140700000 

140800000 

140500000 

141000000 

141100000 

141200000 

141300000 

141400000 

141500000 

141600000 

141700000 

141800000 

141900000 

142000000 

142100000 

142200000 

142300000 

142400000 

142500000 

142600000 

142700000 

142800000 

142900000 

143000000 

143100000 

143200000 

143300000 

143400000 

143500000 

143600000 

143700000 

143800000 



The following list contains the input variables with the dimensions used in the pro- 
gram „ The size of an array is limited to the dimensions stated. The maximum number 
of stations in the x-direction is 20 and the maximum number of stations in the y-direction 
is 10 . 


FORTRAN 

variable 


AEXP 


Symbol 




ALCTAB(IO) a c 


ALP HAT (10) a 

ALSTAB(IO) a- 

D 

ASEXP A s 

BETA (3 


BEXP B 

BSEXP B g 

CE C e 

CKETATB(50) 
CKXITAB(50) k| 

CORDSY 


Description 

Coefficient of the exponential term when the 
Arrhenius expression is used for calculating 
MCDOT 

Aerodynamic -blocking coefficient for heat and 
mass transfer associated with MCDOT, a func- 
tion of time (TTALC) 

Absorptance of surface, a function of temperature 
(TALPHA) 

Aerodynamic -blocking coefficient for heat and 
mass transfer associated with MSDOT, a func- 
tion of time (TTALS) 

Coefficient in the expression for calculating MSDOT 

Determines whether ablation or transpiration 
theory will be used for effect of mass transfer 
on heat transfer; for ablation theory, BETA=1; 
for transpiration theory, BETA=0 

Power of the exponential term in the Arrhenius 
expression for calculating MCDOT 

Power of the exponential term in the expression 
for calculating MSDOT 

Oxygen concentration, by mass, at edge of bound- 
ary layer 

Thermal conductivity in 7] -direction, a function of 
77 (ETATAB) and temperature (TTCKETA) 

Thermal conductivity in £ -direction, a function of 
£ (XITAB) and temperature (TTCKXI) 

Trigger to indicate coordinate system; if curvi- 
linear coordinates, CORDSY=0; if Cartesian 
coordinates, CORDSY=l 
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FORTRAN 

variable 

Symbol 

Description 

CPDP 

e " 
C P 

Specific heat of layer along y=0 

CPP 

c T 
c p 

Specific heat of layer along x=L 

CPTAB(IO) 

°p 

Specific heat, a function of temperature (TTABCP) 

DELTAO(20) 

6 

Initial material thickness, must have L values 

DELTAU 

At 

Initial computing time interval 

DELTMIN 


Minimum value allowed for DELTA 

DTMAX 


Maximum DELTAU which can be used; if no value 
is given, DTMAX=2.0 

ELAMTB(28) 

\ 

Ratio of mass of material removed per unit mass 
of oxygen that reaches the surface, a function of 
pressure (PELAM) and temperature (TTELAM) 

ENDTAU 


Time at which calculation stops 

EPSONE 

e 

Emittance of front surface 

EPSONEP 

e’ 

Emittance of layer along x=L 

EPSONPP 

€" 

Emittance of layer along y=0 

ERRORT 


Acceptable relative error in temperature 

ETATAB(5) 

V 

ETA table for CKETATB 

GAMBAR 


Mean ratio of specific heats behind bow shock 
wave, used only in computation of heating -rate 
distribution around body 

GAMINF 


Ratio of specific heats in free stream, used only 
in computing heating -rate distribution around 
body 

HCOMBTB(28) 

AH c 

Heat of combustion, a function of pressure 
(PHCOMB) and temperature (TTHCOMB) 

HCTAB(28) 

AH S 

Heat of sublimation, a function of pressure (PHC) 
and temperature (TTABHC) 

HETAB(IO) 

H e 

Total free -stream enthalpy, a function of time 
(TTABHE) 

HWTAB(15) 

H w 

Enthalpy of gas at the wall temperature, a function 
of temperature (TTABHW) 
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FORTRAN 

variable 


Symbol 


Description 


IADJUST 


IP LOT 


L 

MACHNO 

MALPHA 

MALPHC 

MALPHS 

MAXITT 


MCP 

MDMAX 

MDOTO(20) 

MHE 

MHW 

MPSTAG 

MQC 

MQR 


Trigger for adjusting heating-rate and pressure 
distributions to shape change; if XADJUST=0, 
QRAT and PRAT are not adjusted; if 
IADJUST/0, QRAT and PRAT will be adjusted 
to shape change 

Trigger for plotting routine; if IPLOT=0, no plots; 
if IPLOT^O, the following plots will be made: 
RSS versus ZS at times indicated in PLTIME 
table; MDOT versus x at each PRFREQ time; 
and T(M,N) versus x indicated in NTP array 
at each PREREQ 

Number of stations in the x-direction 

Free-stream Mach number 

Order of interpolation for ALPHAT 

Order of interpolation for ALCTAB 

Order of interpolation for ALSTAB 

Maximum iteration count; when iteration count 
exceeds this number, DELTAU will be halved 
until DELTAU is less than 1.0E-6, then the 
program will stop and a message will be 
printed 

Order of interpolation for CPTAB 

Maximum expected MDOT; this must be given to 
get a reasonable scale for plots; not needed if 
IPLOT=0 

m Initial mass loss rate at surface, must have 

L values 

Order of interpolation for HETAB 
Order of interpolation for HWTAB 
Order of interpolation for PSTAGTB 
Order of interpolation for QCTAB 
Order of interpolation for QRTAB 
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FORTRAN 

variable 


Symbol 


Description 


MTB 

MW02 

MWSTR 

NALPHA 

NALPHC 

NALPHS 

NCKETA 

NCKXI 

NCP 

NELAM 

NETA 

NHC 

NHCOMB 

NHE 

NHW 

NPELAM 

NPHC 

NPHCOMB 

NPSTAG 

NQC 

NQR 

NTB 

NTP(7) 



M w 


Order of interpolation for TBTAB 

Molecular weight of diatomic oxygen used in 
oxidation equation 

Molecular weight of free stream used in oxidation 
equation 

Number of entries in ALPHAT 

Number of entries in ALCTAB 

Number of entries in ALSTAB 

Number of entries in CKETATB 

Number of entries in CKXITAB 

Number of entries in CPTAB 

Number of entries in ELAMTB 

Number of entries in ETATAB 

Number of entries in HCTAB 

Number of entries in HCOMBTB 

Number of entries in HETAB 

Number of entries in HWTAB 

Number of entries in PELAM 

Number of entries in PHC 

Number of entries in PHCOMB 

Number of entries in PSTAGTB 

Number of entries in QCTAB 

Number of entries in QRTAB 

Number of entries in TBTAB 

Array of seven values which specify the tempera- 
tures to be plotted; NTP(l) = the number of 
temperature rows to be plotted (may be six or 
less); NTP(2) through NTP(7), the row number 
of the temperatures to be plotted. For example, 
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FORTRAN 

variable 


Symbol 


Description 


NXI 

PELAM(4) 

PHC(4) 

PHCOMB(4) 

PLTIME(15) 

PRAT(20) 

PRFREQ 

PSEXP 

PSTAGTB(IO) 

PTMAX 

PTMIN 

QCTAB(IO) 

QRAT(20) 

QRRAT(20) 

QRTAB(IO) 

R(20) 

RIEXP 

RNSI 


P 




q r 

R 


r 

R stag 


NTP(1)=3, NTP(2)=1, NTP(3)=5, NTP(4)=10, 
specifies that three (3) rows of temperature 
will be plotted and these rows are 1, 5, and 10 

Number of entries in XITAB 

Pressure table for ELAMTB 

Pressure table for HCTAB 

Pressure table for HCOMBTB 

Times at which RSS versus ZS, that is, the body 
shape, will be plotted; not needed if IPLOT=0 

Initial ratio of local to stagnation pressure, must 
have L values, not needed if IADJUST^O 

Printing time frequency for output data 

Exponent of pressure term in sublimation 
equation 

Stagnation pressure, a function of time (TTPSTAG) 

Maximum expected value of T, used to get rea- 
sonable scale in plotting, not needed if IPLOT=0 

Minimum expected value of T, used to get rea- 
sonable scale in plotting, not needed if IPLOT=0 

Cold -wall convective heating rate, a function of 
time (TTABQC) 

Initial convective heating -rate distribution must 
have L values, not needed if IADJUST^O 

Radiant heating-rate distribution over body, must 
have L values 

Radiant heating -rate tables, a function of time 
(TTABQR) 

Radius of curvature of base curve at node points, 
must have L values 

Exponent of nose -radius term in MSDOT equation 

Initial nose radius 
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FORTRAN 

variable 


Symbol 


Description 


RO 

RODP 

ROP 

RS(20) 


RSSMAX 


S 

STEBOL 

T(10,20) 

TALPHA(IO) 

TAUO 

TBTAB(IO) 

TDPRIME 

THETA(20) 

TMIN 


P 

P” 

P' 



a 


T 


0 


T PRIME 

TTABCP(IO) 

TTABHC(IO) 

TTABHE(IO) 

TTABHW(15) 

TTABQC(IO) 


Material density 
Density of layer along y=0 
Density of layer along x=L 

Cylindrical radius from body axis of symmetry to 
node points on the base curve, must have 
L values 

Maximum expected value of RSS, used to get a 
reasonable scale for plots, not needed if 
IP LOT =0 

Number of stations in y-direction 
Stefan -Boltzmann constant for radiation 
Initial temperature, must have S*L values 
Temperature table for ALPHAT 
Initial time 

Temperature to which back surface is radiating, 
a function of time (TTABTB) 

Thickness of layer along y=0 

Angle (in degrees) less than or equal to 90° 
between RS and R, must have L values 

Minimum temperature value; if TMIN^O and a 
computed temperature goes below TMIN, the 
temperature will be set equal to TMIN; if 
TMIN=0, no restraint will be made on the com- 
puted temperatures 

Thickness of layer along x=L 

Temperature table for CPTAB 

Temperature table for HCTAB 

Time table for HETAB 

Temperature table for HWTAB 

Time table for QCTAB 
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FORTRAN 

variable 

Symbol 

Description 

TTABQR(IO) 


Time table for QRTAB 

TTABTB(IO) 


Time table for TBTAB 

TTALC(IO) 


Time table for ALCTAB 

TTALS(IO) 


Time table for ALSTAB 

TTCKETA(IO) 


Temperature table for CKETATB 

TTCKXI(IO) 


Temperature table for CKXITAB 

TTELAM(7) 


Temperature table for ELAMTB 

TTHCOMB(7) 


Temperature table for HCOMBTB 

TTPSTAG(IO) 


Time table for PSTAGTB 

XITAB(5) 

1 

Table of values of CKXITAB 

XO 

x b 

Length of base curve 

XORDER 


Order of oxidation 

ZS(20) 


Initial distance from the initial stagnation point 
to RSS along body axis of symmetry, must 
have L values 

ZSMAX 


Maximum expected value of ZS, used to get rea 
sonable scale for plotting RSS versus ZS, not 
needed if IPLOT=0 < 



Output 


Examples of output data are given in appendix B. The input data are printed at the 
beginning of the output listing in the same order in which they appear in the NAMELIST 
statement. Then the identification card is printed. Headings and interpretations are as 
follows: 

Heading Description 

TAU Time at which the calculations were made 

DELTAU The computing time interval 

QC Convective heating rate 

QR Radiant heating rate 
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Heading 


HE 

T(S,1) 


TEMPERATURE 

(M,N) 


ETA 

X 

MDOT(N) 

MCDOT(N) 

DELTA(N) 

QRAT(N) 

PRAT(N) 

QS(N) 

QRR(N) 

QCOMB(N) 

NO.ITER.COL. 

NG.ITER.ROW. 

TOTAL NO.ITER. 

IROCOL 


Total free -stream enthalpy 

Temperature at time r - At; this value can indicate whether 
a reasonable At is being used; by observing this value 
and the value at t, unusual behavior might indicate the 
need for a smaller At. 

Temperatures; to locate the station read ETA to the left and 
x above the temperature column; up to 15 temperatures are 
printed on one line; if more columns have been used, the 
remaining temperatures will be printed on the next line 

Dimensionless y values, printed in the first column on the 
left side of the page 

Length along base curve from stagnation point to the station, 
printed in the second column and reading from left to right 

Surface mass loss rate at station n 

Mass loss due to oxidation at station n 

Material thickness at station n 

Ratio of local heating to stagnation heating at station n 
Ratio of local pressure to stagnation pressure at station n 
Net heat input at station n 
Reradiation at station n 

Heat due to combustion for oxidation at station n 

Number of iterations for the previous column solution 

Number of iterations for the previous row solution 

Total number of iterations from the beginning of the problem 

Tells at which solution the printout was made; value of 1 indi- 
cates column solution; 2, row solution 


Diagnostics 

The program has several automatic stops to avoid the waste of computer time on 
problems which appear to be having computational difficulties. These stops are 

(1) DELTA < DELTMIN: If any thickness DELTA becomes less than the input 
DELTMIN a normal printout is made and the program will stop. 
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(2) Negative temperature: If any temperature becomes negative, a normal printout 
is made and the program will stop. 

(3) DELTAU < 1.0E-6: If the computing time interval DELTAU becomes less than 
1.0E-6, the message TEMPERATURE ITERATION DOES NOT CONVERGE will be printed. 
The current estimated temperatures are printed, a normal printout is made, and the pro- 
gram will stop. 

(4) Iteration count exceeded: If the maximum iteration count input MAXITT is 
exceeded and the calculation is a row solution, the computing interval cannot be halved. 
The message THIS IS A ROW SOLUTION, DELTAU CANNOT CHANGE is printed. The 
current estimated temperatures are printed, a normal printout is made, and the program 
will stop. 

(5) Temperature diverging: If any temperature begins diverging, the message 

TEMPERATURE IS DIVERGING WHY is printed. The current estimated temper- 

atures are printed, a normal printout is made, and the program will stop. 

Whenever these diagnostics appear, the input should be checked to make sure that 
all initial conditions have been given. Check all input tables for any discontinuities. 
Negative temperatures may result from oscillations caused by time intervals which are 
too large. High values of MDOT and rapid changes of heat input with time may require 
smaller time intervals for computational purposes. 

SAMPLE CASES 

Three sample cases are presented to illustrate the operation of the computer pro- 
gram. All the cases are for ablating bodies of different geometries: a hemisphere, a 
hemispherically blunted cone, and a right -circular cylinder. A listing of the input data 
and a sample of the output data for each case are shown in appendix B. 

Computer-generated curves of some of the output from the sample cases are shown 
in figures 1, 2, and 3. The curves show body shape change due to ablation, histories of 
mass -transfer rate over the surface of the bodies, and selected temperature histories. 
The body shape is plotted at each time listed in the input PLTIME. The mass-transfer 
rates over the surface and the temperatures along the rows specified by the input NTP 
are plotted at each printing frequency for the output data. 
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The computing time depends on the accuracy desired; the boundary condition, that 
is, the heating -rate history; and the number of node points. The computational times for 
the sample cases are 136 seconds for the hemisphere, 312 seconds for the right -circular - 
cylinder, and 150 seconds for the hemispherically blunted cone. These cases have not 
been optimized with respect to time and, therefore, may run in shorter periods of time. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., September 3, 1971. 
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APPENDIX A 


LANGLEY LIBRARY SUBROUTINES 
Subroutine FTLUP 


Language : FORTRAN 

Purpose: Computes y = F(x) from a table of values using first- or second-order interpolation. 

An option to give y a constant value for any x is also provided. 

Use : CALL FTLUP(X, Y, M, N, VARI, VARD) 

X The name of the independent variable x„ 

Y The name of the dependent variable y = F(x). 

M The order of interpolation (an integer) 

M = 0 for y a constant. VARD(I) corresponds to VARI(I) for 

1=1,2,. , ., N. For M = 0 or N 2 1, y = F(VARI(1)) for any value of x. 

The program extrapolates. 

M = 1 or 2. First or second order if VARI is strictly increasing (not equal). 

M = -1 or -2. First or second order if VARI is strictly decreasing (not equal). 

N The number of points in the table (an integer) . 

VARI The name of a one-dimensional array which contains the N values of the independent variable. 

VARD The name of a one-dimensional array which contains the N values of the dependent variable. 

Restrictions : All the numbers must be floating point. The values of the independent variable x in the 
table must be strictly increasing or strictly decreasing. The following arrays must be dimensioned by 
the calling program as indicated: VARI(N), VARD(N). 

Accuracy : A function of the order of interpolation used. 

References : (a) Nielsen, Kaj L.: Methods in Numerical Analysis. The Macmillan Co., c.1956, pp. 87-91. 

(b) Milne, William Edmund: Numerical Calculus. Princeton Univ. Press, c.1949, pp. 69-73. 

Storage : 430g locations. 

Error condition: If the VARI values are not in order, the subroutine will print TABLE BELOW OUT OF 
ORDER FOR FTLUP AT POSITION xxx TABLE IS STORED IN LOCATION xxxxxx (absolute). It then 
prints the contents of VARI and VARD, and STOPS the program. 

Subroutine date: September 12, 1969. 
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APPENDIX A - Continued 


Subroutine DISCOT 

Language : FORTRAN 

Purpose : DISCOT performs single or double interpolation for continuous or discontinuous functions. 
Given a table of some function y with two independent variables, x and z, this subroutine performs 
K x th- and K z th-order interpolation to calculate the dependent variable. In this subroutine all single- 
line functions are read in as two separate arrays and all multiline functions are read in as three 
separate arrays; that is, 


x i 

(i= 1,2, . . 

-,L) 

y j 

(3 = 1,2,. . 

M) 

44 

N 

CM 

T— 1 
11 

44 

-,N) 


Use : CALL DISCOT (XA, ZA, TABX, TABY, TABZ, NC, NY, NZ, ANS) 

XA The x argument 

ZA The z argument (may be the same name as x on single lines) 

TABX A one-dimensional array of x values 

TABY A one -dimensional array of y values 

TABZ A one-dimensional array of z values 

NC A control word that consists of a sign (+ or -) and three digits. The control word is formed 

as follows: 

(1) If NX = NY, the sign is negative. If NX * NY, then NX is computed by DISCOT as 
NX = NY/Nz, and the sign is positive and may be omitted if desired. 

(2) A one in the hundreds position of the word indicates that no extrapolation occurs above 
z max . w Wh a zero in this position, extrapolation occurs when z > z max . The zero 
may be omitted if desired. 

(3) A digit (1 to 7) in the tens position of the word indicates the order of interpolation in 
the x-direction. 

(4) A digit (1 to 7) in the units position of the word indicates the order of interpolation in 
the z-direction. 

NY The number of points in y array 

NZ The number of points in z array 

ANS The dependent variable y 
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APPENDIX A - Continued 


The following programs will illustrate various ways to use DISCOT: 


CASE I: Given y = f(x) 

NY = 50 

NX (number of points in x array) = NY 
Extrapolation when z > z max 
Second-order interpolation in x-direction 
No interpolation in z-direction 
Control word = -020 
DIMENSION TABX (50), TABY (50) 

1 FORMAT (8E 9.5) 

READ (5,1) TABX, TABY 
READ (5,1) XA 

CALL DISCOT (XA, XA, TABX, TABY, TABY, -020, 50, 0, ANS) 

CASE II: Given y = f(x,z) 

NY = 800 
NZ = 10 

NX = NY/NZ (computed by DISCOT) 

Extrapolation when z > z max 
Linear interpolation in x-direction 
Linear interpolation in z-direction 
Control word = 11 

DIMENSION TABX (800), TABY (800), TABZ (10) 

1 FORMAT (8E 9.5) 

READ (5,1) TABX, TABY, TABZ 
READ (5,1) XA, ZA 

CALL DISCOT (XA, ZA, TABX, TABY, TABZ, 11, 800, 10, ANS) 

CASE III: Given y = f(x,z) 

NY = 800 
NZ = 10 
NX = NY 

Extrapolation when z > z max 
Seventh-order interpolation in x-direction 
Third-order interpolation in z-direction 
Control word = -73 

DIMENSION TABX (800), TABY (800), TABZ (10) 

1 FORMAT (8E 9.5) 

READ (5,1) TABX, TABY, TABZ 
READ (5,1) XA, ZA 

CALL DISCOT (XA, ZA, TABX, TABY, TABZ, -73, 800, 10, ANS) 

CASE IV: Same as Case IE with no extrapolation above z max . Control word = 

CALL DISCOT (XA, ZA, TABX, TABY, TABZ, -173, 800, 10, ANS) 


-173 
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APPENDIX A - Continued 


Restrictions : See rule (5c) of section "Method" for restrictions on tabulating arrays and discontinuous 
functions. The order of interpolation in the x- and z-directions may be from 1 to 7. The following 
subprograms are used by DISCOT: UNS, DISSER, LAGRAN. 


Method : Lagrange's interpolation formula is used in both the x- and z-directions for interpolation. This 
method is explained in detail in reference (a) of this subroutine. For a search in either the x- or 
z-direction, the following rules are observed: 

(1) If x < Xj, the routine chooses the following points for extrapolation: 

x i ,x 2> • . x k+1 and y 1 >y 2 > • • •> y^+i 

(2) If x > x n , the routine chooses the following points for extrapolation: 

x n-k ,x n-k+l’ ' * - x n and Wa-k+1" ’ y n 

(3) If x S x n , the routine chooses the following points for interpolation: 

When k is odd, 


2 2 


4 k+1 and y k+l’ y k+1 ’ 

i-iSii+k i-£E±i i-S+i+1 

2 2 2 


• ? y v 1 1 

i-£±i + k 

2 


WHhen k is even. 


‘-r-K ■ 


k and y k ,y k 


i-t +k 


i-£ i-5+1 
2 2 


y k 

i -i +k 


(4) If any of the subscripts in rule (3) become negative or greater than n (number of 
points), rules (1) and (2) apply. When discontinuous functions are tabulated, the indepen- 
dent variable at the point of discontinuity is repeated. 

(5) The subroutine will automatically examine the points selected before interpolation and if 

there is a discontinuity, the following rules apply. Let x^ and be the point of 

discontinuity. 

(a) If x S x^, points previously chosen are modified for interpolation as shown: 

x d-k’ x d-k+l’ • • -> x d y d-k’ y d-k+l> • • -> y d 

(b) If x > x^, points previously chosen are modified for interpolation as shown: 

x d+l ,x d+2’ ' • *’ x d+k and y d+l’ y d+2’ y d+k 

(c) When tabulating discontinuous functions, there must always be k+1 points above 
and below the discontinuity in order to get proper interpolation. 

(6) When tabulating arrays for this subroutine, both independent variables must be in 
ascending order. 
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APPENDIX A - Concluded 


(7) In some engineering programs with many tables, it is quite desirable to read in one 
array of x values that could be used for all lines of a multiline function or different 
functions. Even though this situation is not always applicable, the subroutine has been 
written to handle it. This procedure not only saves much time in preparing tabular 
data, but also can save many locations previously used when every y coordinate had 
to have a corresponding x coordinate. Another additional feature that may be useful 
is the possibility of a multiline function with no extrapolation above the top line. 

Accuracy : A function of the order of interpolation used. 

Reference : (a) Nielsen, Kaj L.: Methods in Numerical Analysis. The Macmillan Co., c.1956. 

Storage : 555g locations. 

Subprograms used : UNS 40g locations. 

DISSER 110g locations. 

LAGRAN 55g locations. 

Subroutine date : August 1, 1968. 
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The sample output listing for a teflon hemisphere is given below 
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The sample input listing for a graphite hemisphere -cone is given below: 
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The sample output listing for a graphite hemisphere-cone is given below 
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The sample input listing for a right -circular cylinder is given below 
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MDMAX=.0?=;, PTMaX= 850. PTMIN = 350* RSSMAX=.l. ZSM AX= • 05 . 
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ITER. COL.= 1 NO. ITER. ROW= 1 TOTAL NO. ITER.= 17A IROCOL= 
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Figure 3.- Computer-generated profile, mass loss, and temperature 
histories for a right-circular cylinder. 
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